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ABSTRACT 

We present a semi- analytic model for Cold Dark Matter halo substructure that can be used as a 
framework for studying the physics of galaxy formation and as an ingredient in halo models of galaxy 
clustering. The model has the following main ingredients: (1) extended Press-Schechter mass accretion 
histories; (2) host halo density profiles computed according to the trends observed in cosmological simu- 
lations; (3) distributions of initial orbital parameters of accreting subhalos measured in a high-resolution 
simulation of three Milky Way-size halos; and (4) integration of the orbital evolution of subhalos includ- 
ing the effects of dynamical friction and tidal mass loss. We perform a comprehensive comparison of the 
model calculations to the results of a suite of high-resolution cosmological simulations. The comparisons 
show that subhalo statistics, such as the velocity and mass functions, the radial distributions, and the 
halo occupation distributions agree well over three orders of magnitude in host halo mass and at various 
redshifts. We find that both in the simulations and in our model the radial distributions of subhalos 
are significantly shallower than that of the dark matter density. The abundance of subhalos in a host 
is set by competition between tidal disruption and new accretion. Halos of high mass and halos at high 
redshift tend to host more subhalos because the subhalos have, on average, been accreted more recently. 
Similarly, at a fixed mass and epoch, halos that formed more recently host a larger number of subhalos. 
Observed "fossil groups" may represent an extreme tail of this correlation. We find a related correlation 
between host halo concentration and satellite abundance at fixed host mass, A'^gat oc c~°:, where a changes 
with redshift and host-to-subhalo mass ratio. Lastly, we use our substructure model to populate host 
halos in one of the high-resolution cosmological simulations, replacing the actual subhalos resolved in this 
simulation and using the host mass as the only input for the model calculation. We show that the result- 
ing correlation function of such a hybrid halo ensemble is indistinguishable from that measured directly 
in the simulation. This supports one of the key tenets of the standard halo model — the assumption 
that the halo occupation distribution is statistically independent of the host halo environment. 

Subject headings: cosmology: theory, galaxies: formation, large-scale structure of universe 



1. INTRODUCTION 

Numerical simulations of structure formation set 
within the co l d dar k matter (CDM) p aradig m (see 
IWhite k Reei Il978t iBlumenthal et all Il984li have 
revealed that virialized, dark matter halos are rife 
with disti nct, 
los 
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2004t iDiernand. Moore. fc Stadell 12004 iGao et al.ll2004al: 
Reed et al.ll2004^ . It is tempting to associate subhalos. 



particula rly large subhalos, w ith gal axies in group s and 
clusters (iKlvpin et alj Il999al: Krav tsov fc KlvpinI [19991 
IColm et aL 1999t Moore et al.i il999l). Le nding support 
to this conjecture, Kravtsov et al.l (|2004l ) used natural 
assumptions to connect galaxy luminosity to subhalos 
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and showed that the clustering properties of subhalos and 
their host halos mat ch the observed cluste ring properties 
of galaxies (see also iNevrinck et al.l |2004|) . This implies 
that the occupation distribution of subhalos must be 
quite similar to th e occupation distribution of galaxies 
in host halos (e.g. iBerlind fc WeinberdI l2002t ) and that 
luminous satellite galaxies are likely associated with 
massive dark matter subhalos. A corollary of this is that 
if one understands the physical processes that govern 
subhalo occupation, one gains significant insight into the 
physics of galaxy clustering. 

The abundances of subhalos within host dark matter 
halos are determined through the interplay of several im- 
portant and relatively distinct physical processes. If one 
imagines, for simplicity, the development of a subhalo pop- 
ulation as a sequence, the first process to consider is the 
merging of halos to form ever larger systems. In the pre- 
vailing CDM paradigm of hierarchical structure growth, 
halo merger rates are governed by the initial density field 
and the global cosmological model, and tend to be more 
rapid at early times. The second step is to consider the 
evolution of merged subhalos in the potential of the host 
halo. In contrast to merger statistics, subhalo evolution 
is influenced by local processes in the high-density envi- 
ronment of the host, like dynamical friction, tidal mass 
loss, and heating by violent interactions with other struc- 
tures. These processes reduce the amount of distinct sub- 
structure within host systems over timescales of order the 
local dynamical time. At early times, when the merger 
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timescale is short compared to the dynamical time, we ex- 
pect high substructure counts within individual host ha- 
los. At late times, we expect substructure to be reduced, 
as subhalos are packed into more and more dense environ- 
ments inside larger and larger host systems. All of these 
physical processes are purely gravitational, but have a sig- 
nificant effect on the small-scale clustering of subhalos. 
By implication, galaxy clustering statistics are governed 
by these processes as well. 

This paper is the first of a series of studies of the physical 
processes that govern galaxy clustering statistics. In the 
context of this series, this paper represents a description of 
our theoretical methodology, but the content of this paper 
is considerably more general. We present a semi-analytic 
model that can be used to make predictions for the sub- 
structure populations of dark matter halos. The model be- 
gins with the computation of the merger histories of dark 
matter halos and approximates the subsequent effects of 
dynamical friction and tidal mass after accretion. In the 
regime where they are commensurable, we make detailed 
comparisons of this model to the results of high-resolution 
iV-body simulations to assess the accuracy and applicabil- 
ity of the model (see §[!]). We emphasize that the model we 
present in this first paper includes only the effect of gravity 
and neglects any effects of baryon cooling and condensa- 
tion. Our strategy in so doing is to build upon a tractable 
model that can be robustly be tested with A^-body simu- 
lations in a variety of ways. In this way, the importance 
of additional processes and approximations can be ascer- 
tained in a well-controlled way and we will pursue these 
aims in the subsequent papers of this series. 

There are advantages to using an analytic model that 
has been tested against existing simulations in a wide 
range of applications. First, an analytic model can quickly 
generate a statistically-large number of realizations of sub- 
halo populations within host halos of various masses with 
a relatively small computational effort. As a result, such 
a model can complement direct simulation where limited 
dynamic range severely restricts the size of any host halo 
sample in a study of ha l o subs tructure. Along these hues, 
llslam. Tavlor. and Silkl (|2003l ) used a semi-analytic ap- 
proach to study a population of black holes, which may ex- 
ist as remnants of the first stars, in the Milky Way (MW) 
halo. This study would have required a dynamic range 
that is unachievable even with dissipationless numerical 
simulations. In addition, such a model can be used to 
study the effects of non-standard cosmological inputs on 
subhalo population s over a wide r ange of parameter space. 
IZentner fc Bullockl (|2003l . IZB03I hereafter) used such an 
approach to estimate the influence of modifications to the 
primordial power spectrum on the populations of satellite 
halos and the prospects for using satellite halo populations 
to constrain cosmology. 

More closely tied to the specific aims of the papers in 
this series, analytic substructure models can provide an 
important ingredient for modeling the large-scale clus- 
tering of galaxies by predicting halo occupation distri- 
butions. Coupled with a relatively low-resolution, large- 
volume simulation, such a model can be used to calcu- 
late galaxy clustering statistics, including environmental 
dependences, in a way that explicitly accounts for com- 
plicated non-linear effects such as the bias of host halos 
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ll998HSheth. Mo. fc Tori^ll999l : ISeliak fc Warrenll2004D 
Alternatively, a subhalo model can be utilized in conjunc- 
tion with an analytic h alo model (see, e.g. Scliak 2 0001 : 
IScoccimarro et alll2001[ ) to make analytic predictions for 
clustering statistics rapidly, for a wide range of the space 
of cosmological parameters and for various assumptions 
about galaxy formation. We demonstrate this utilization 
explicitly in § 14.61 Importantly, when used in these ways, 
our subhalo model provides a framework for rapidly testing 
hypotheses about the physics of galaxy formation against 
the observed clustering properties of galaxies. Of equal 
importance, the results of the semi-analytic model predic- 
tions are easy to dissect and understand in a physical way. 
We explore these avenues in the subsequent papers of this 
series. 

The outline of this paper is as follows. In § [2J we 
briefly describe three sets of A^-body simulations that 
we use to test our analytic subhalo model over a broad 
range of host halo properties. In § [31 we describe all of 
the ingredients of the substructure model. The model 
is an extension of the work of ZBOSi and we pay par- 
ticular attention to differences and improvements upon 
that work. As we discuss below, our mo del is similar in 
many respect s to the model developed by iTavlor fc Babull 
()2001l l2004a[). and shares common ingredien ts wit h the 
models of iBenson et al.l (I2002D. [Ogliri fc Le a (|2004D . and 
Ivan den Bosch. Tormen. fc Gioccolir ( 2004bl ). We present 
our results and compare them to the results of high- 
resolution numerical simulations in § |4l In this section, we 
also discuss our results in the context of many recent obser- 
vational and theoretical developments and more generally 
describe the utility of our substructure model. We draw 
conclusions in § [5] In a forthcoming companion paper (A. 
A. Berlind et al., in preparation), we use this model to ad- 
dress the origin of the nearly power-law galaxy correlation 
function and several other features of galaxy clustering 
statistics. 

Throughout this paper we refer to self-bound subunits 
that are within the virial radius of another, larger halo 
as subhalos or satellites. We refer to halos that are not 
contained within another, larger halo, and therefore are 
not classified as subhalos, as host halos. We use the term 
halo to refer to both host halos and satellite halos. 



2. NUMERICAL SIMULATIONS 

The logic behind our use of numerical simulations is as 
follows. By their nature, semi-analytic models attempt to 
model complicated, non-linear dynamical processes using 
simple prescriptions and rely on inputs distilled from nu- 
merical simulations. Some of the inputs for the evolution 
of subhalos is provid ed by controlled (non-cosmo l ogical) 
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20031 : Kazantzidis et al. l2004bD . However, cosmological 
simulations are still needed to specify the initial orbital 
parameters of accreted subhalos and the structure of 
cosmological halos. We assume halo density profiles 
derived from cosmological simulations for a wide range 
of masses and cosmologies as described below. We 
determine the initial orbital parameters of subhalos from 
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a high-resolution simulation of three nearly MW-size 
hales. We subsequently compare the results of our 
model with the subhalos in simulations of a very high 
resolution cluster-size halo and a large number of halos in 
a cosmological simulation with uniform mass resolution 
and a different power spectrum normalization. These 
comparisons over a wide range of halo masses and different 
power spectrum normalizations are thus non-trivial tests 
of the applicability of our model. 

In both our simulations and our semi-analytic model- 
ing, we quantify the size of host dark matter halos by 
their virial masses M^ir or, equivalently, their virial radii 
-Rvir- We define the halo virial radius as the radius within 
which the mean density is equal to the virial overdensity 
Avir, multiplied by the mean matter density of the uni- 
verse pM, so that Mvir = 47rpMAviri?vir/3- For the sim- 
ulated halos, this radius is centered on the particle with 
the highest local density. The quantity Avir can be esti- 
mated from the sphe rical top-hat collapse approximation 
(|Eke. Navarro, fc EY enk 1998). We compute Ay j, - usin g the 
fitting function provided bv iBrvan fc NormanI (|1998l ). In 
the ACDM cosmology that we consider, Avir (2: = 0) ~ 337 
and tends toward the standard CDM (i.e. ^Im = 1) value 
Avir 178 at high redshift {z> 1). 

All of our simulations were performed with the 
Adaptive Refinement Tree (ART) iV-body code (see 



iKravtsovT Klvpin. fc Khokhlovlll997l : lKravtsovll99l . The 

simulated Galax y-size halos are th e three halos dis- 
cussed in de tail in'Klvpin et al.' (2001', hereafter 'KOT) and 
iKravtsov. Gn cdin. & Klypin (2004, hereafter KGKOl). 
Briefly, iKOll used the ART code to model the evolution 
of three MW-size halos in a standard "concordance," flat, 
ACDM cosmology with = 1 — ^^A = 0.3, h = 0.7, and 
as = 0.9 in a comoving box of size 25 h^^Mpc on a side. 
The simulation began with a uniform 256'^ grid covering 
the entire computational volume. Higher force resolution 
was achieved about collapsing structures by recursive re- 
finement of all such regions using an adaptive refinement 
algorithm. The grid cells were refined if the particle mass 
contained within them exceeded a certain specified thresh- 
old value. 

The multiple mass resolution technique was used to 
set up the initial conditions. A Lagrangian region cor- 
responding to two virial radii about each halo was re- 
sampled with the highest resolution particles of mass 
nip — 1.2 X 10^ H^^Mq, corresponding to 1024'^ particles 
in the box, at the initial timestep, Zi = 50. The high mass 
resolution region was surrounded by layers of particles of 
increasing mass with a total of five particle species. Only 
the regions containing the highest-resolution particles were 
adaptively refined. The maximum level of refinement in 
the simulations corresponded to a peak formal spatial res- 
olution of approximately 100/i~^pc. The three MW -size 
halos, referred to as Gi, G2, and G3 in |KGK04 have 
virial masses of 1.13 x 10^^ /^-iMq, 1.14 x lO^^ /^-iMq, 
and 1.45 x 10^^ H'^Mq respectively, and each of them is 
resolved with approximately ~ 10^ highest-resolution par- 
ti cles w ithi n their vi rial radii. Further details can be found 
in iKOlland lKGKOi . 

In addition to these three MW-size halos, we com- 
pare our model results to the results of a high-resolution 
simulation of a cluster-size halo of virial mass ^ 1.7 x 



10^"* /i~^Mq. The cluster was simulated in a comoving 
box of 80 /i^^Mpc on a side using the same multiple mass 
resolution technique described above in the same concor- 
dance ACDM cosmology, as described in Tasitsio mi et al.l 
(|2004l ). The peak formal resolution was 0.6 h~^kpc and 
the mass of the highest resolution particles was rup = 
3.95 X 10*" H^^Mq. The cluster was resolved with f« 
4.4 X 10^ particles within its virial radius. This high- 
resolut ion simulated cluster h alo was referred to as CL2 
HR bv lTasitsiomi et al.l (|2004f) . 

Lastly, we compare our model predictions to the sub- 
halo populations in a cosmological simulation with uni- 
form mass resolution containing a large number of host 
halos over a wide range masses. This simulation fol- 
lowed the evolution of 512^^ particles in an 80 h^^Mpc 
comoving box, corresponding to a particle mass of nip = 
3.16 X 10^ /i~^Mq. The cosmological model was the same 
ACDM cosmology as above, but with a "low" power spec- 
trum normalization of CTg = 0.75. The peak formal spa- 
tial resolution of the simulation was 1.2/i~^kpc comoving. 
At 2; = the simulation box contains ~ 1.2 x 10^ host 
halos with Mvir > 10" h'^MQ and « 3 x 10^ subhalos 
with a bound mass > 10^^ h~^MQ. We use this simu- 
lated halo and subhalo population to perform a compre- 
hensive test of our mo del. This simulation was described in 
IKravtsov et al.l (|2004D . where the simulation was referred 
to as ACDMso. 

In all three of the numerical simulations described 
above, halos and subhalos were identified us ing a variant 
of the Bound Density Maxima Algorithm (|Klvpin et al.l 
1999b). The first steps of the halo-finding algorithm are 
to compute the local overdensity at each particle position 
using a smoothing kernel of 24 particles and to find the lo- 
cations of local maxima in the density field. Starting with 
the highest density particle and proceeding toward lower 
density, each density peak is marked as a potential halo 
center and then surrounded by a sphere of fixed radius 
ffind- AH particles within rtind are excluded from further 
consideration as potential halo centers. The search radius 
parameter rfind , is set according to the size of the smallest 
object that we aim to identify. For the simulated MW-size 
halos this was set to r^nd = 10/i~^kpc, while for the high- 
resolution cluster simulation and the ACDMgo simulation 
we set rfind = 50/i~^kpc. After identifying potential halo 
centers, we analyze the surrounding particl es and itera- 
tively remove particles that are unbound fsee lKlvpin et al.l 
Il999bfl . All remaining bound particles are then used to 
compute halo/subhalo properties such as mass M, the cir- 
cular velocity profile Vdrdr) = \/GM{< r)/r, and the 
peak circular velocity Vinax- 

For subhalos, the outer boundary or outer radius of 
the system can be somewhat ambiguous and dependent 
upon a particular definition. Our convention is to adopt 
a truncation radius rt for subhalos, at which the loga- 
rithmic slope of the density profile constructed from the 
bound particles becomes greater than a critical value of 
dln(/9)/dlnr = —0.5. The iterative removal of unbound 
particles is imperfect, and for subhalos some fraction of 
unbound particles with nearly constant density is left on 
the outskirts of the subhalo. Some of these particles are 
from the diffuse mass of the background host halo and 
often some are unbound particles from the tidal debris 
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of the subhalo itself. Generally, the density profiles of 
the tidally-truncated subhalos have an outer slope steeper 
than dln(/9)/dlnr = —3 until the radius where the den- 
sity of background particles becomes comparable to the 
density of the residual background particles not removed 
by the unbinding procedure. At this radius, the profile 
flattens significantly and this is the radius that we iden- 
tify as the subhalo truncation radius. The above criterion 
is thus based on the fact that we do not expect the den- 
sity profiles of CDM halos to be shallower than this and, 
empirically, this definition of truncation radius is approxi- 
mately equivalent to the radius at which the density of the 
bound particles is equal to the background host halo den- 
sity. We note that the outer profiles of subhalos generally 
fall off with dln(p)/dlnr < —3, so the bound mass con- 
verges well before the truncation radius is reached. Con- 
sequently, any uncertainty in the truncation radius docs 
not translate into an appreciable error in the bound mass 
of the subhalo. However, throughout most of this work, 
we quantify the size of satellite halos by their maximum 
circular velocities T^nax, because this quantity is measured 
more robustly and is not subject to the same ambiguity as 
a particular mass definition. This makes our results easier 
to compare to the results of other researchers using differ- 
ent subhalo identification algorithms and mass definitions. 

3. A MODEL FOR HALO SUBSTRUCTURE 

In order to determine the substructure properties of 
a dark matter halo we must model its mass accre- 
tion history as well as the orbital evolution of all ac- 
creted subhalos once they are incorporated into the 
host system. We model substructure using a semi- 
analytic technique that incorporates simplifying approx- 
imations and empirical relations observed in numeri- 
cal simulations. The model is an update d and im- 
pr oved version of the model des c ribed in IZB03I and 
in iKous hiappas. Zentne r. fc Walked ([2 004') and draws on 
the ea rlier work of Bullock. Kravtsov. fc Wein berg (2000, 
doOli). The model of IZBoJ and the improved model that 
we describe here are similar in many respects to the models 
develo pe d byl^ ylor & Babul ([20011, [^04a), Benson et al 
(2002), OgurifcLea (|2004l) . and Tvan den Bosch et al 



2004h ). In this section, we give a brief description of 



the specific model that we use, highlighting additions and 
improvements. We begin with a brief review of the gen- 
eration of merger histories followed by a discussion of our 
assumptions about the density profiles of CDM halos and 
subhalos. We conclude this section with a description of 
our initial conditions for subhalo orbits and our treatment 
of subhalo dynamics within the host potential. 

3.1. Merger Trees 

The first step in determining the properties of halo sub- 
structure is to model the mass accretion history of the host 
dark matter halo. A merger tree that approximates many 
of the results of A^-body simulations can be constructed 
from the linear power spectrum of density perturbations 
using a statistical Monte Carlo technique b ased on the 
extended Press-Schechter (EPS) fo rmalism (iBond et al 
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to form a host halo of a given mass at a given redshift. In 
this way we are able to track both diffuse mass and sub- 
halo accretion via mergers for each host. By generating 
a large number of such merger trees, we can attempt to 
sample the diverse mass accretion histories that result in 
a host halo of a given fixed mass at the time at which we 
are interested in "observing" the halo (most often, z = 0). 

Several authors have explored implementations of spe- 
cific variations of the EPS formalism to construct realistic 
merger trees (e.g. llcgllSheth fc LemsonlHoollCole et all 
20 001). We employ the part icular algorithm advocated by 



Somerville fc KolattI (|1999( ). which conserves mass explic- 
itly. Let a{AI) represent the rms fluctuation amplitude in 
the density field smoothed on a scale R such that the mean 
mass contained within a sphere of radius R is M, that is 
M = AirpyiR^ /"i where pM is the mean mass density of the 
Universe. Following the notation established in LC9^, let 
us denote S{M) = a^{M), AS = S{M) - S{M + AM), 
'w{t) = Sc{t), and Sw = w{t) — w{t + At), where 6c{t) is 
the linear overdensity for collapse at time t in our choic e 
of cosmological model {Sc w 1.68, see lLC93l : lWhitel[l996h . 
With these definitions, the probability that a halo of mass 
M at time t accretes an amount of mass associated with a 
change AS in a timestep corresponding to 5w is given by 



P{AS, Sw)d{AS) 



6w 



-jSw) 

2¥(AS')3/2 2A5 
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19941) . Using this method, we generate a list of the masses 



and accretion redshifts of all subhalos that have merged 



Following rSomerville fc KolattI (|l999t ). we generate merger 
histories by starting at the redshift that we wish to ob- 
serve the final halo (usually z — 0) and step backward 
in time. At each timestep, we select halo progenitors 
by drawing values of AM from the distribution of Equa- 
tion In this way, we generate a list of progeni- 
tor masses and accretion redshifts at each time step. If 
the minimum halo mass that we wish to track is Mmin, 
then in order to reproduce the conditional mass func- 
tion of EPS theory we must choose a timestep such that 
6w = A//tsA^mindS'(A/„ii„)/dM with /ts < 1. In prac- 
tice, we choose fts = 10~^ as an adequate compromise 
between accuracy and computation time. We retain all 
information about mergers with AM > Mmin and treat 
events with AM < Mj^m as diffuse mass accretion. In 
what follows, we employ several different values of Mmin, 
depending upon the minimum halo size that we wish to 
track in each case. 

3.2. Halo Mass Density Profiles 

After accretion, each subhalo is subject to various dy- 
namical processes as it orbits within the host halo. The 
amount of mass that remains bound to the subhalo and 
whether or not the subhalo survives as a distinct, self- 
bound system depend upon the density structure of both 
the subhalo and the host halo. In this subsection we de- 
scribe our assumptions about the mass density profiles of 
CDM halos. 

We characterize the size of host halos by their 
virial masses Afvir and virial radii i?vir, as described 
in § [5] In the interest of simplicity, we model all 
halos with the sp herically-averaged density profile of 
iNavarro. Frenk. fc White (,1997. . .NE^ hereafter): 
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1 



(2) 



where Cvir is the concentration parameter, which describes 
the scale radius = -Rvir/cvin where the logarithmic slope 
of the profile, din p/dlnr = —2. The normalization is 
fixed by the requirement that the mass interior to i?vir 
be equal to Afvir- The NFW profile can be rewritten as 
a circular velocity profile V^i^^{r) = VXcvirg{x)/xg{c^ir), 
where x = r/rg, the function g{y) = ln(l + y) — y/{l + y), 
and Vvir = \/GA'IviT/Rvii- is the halo virial velocity. The 
maximum of the circular velocity profile occurs at a radius 
J-max ^ 2.16rs, at a value of V^^^ ~ 0.22V;?i^Cvir/ff(cvir)- 

The concentrations of dark matter halos observed in 
simulations are tightly correlated with their mass ac cretion 
histo ries (jWcchslcr et al. 2002, hereafter W02; Zhao et al.l 
|2003|) . We account for the correlation between the con- 
centrations of host ha los an d their mass accretion histories 
using the proposal of lW02l . Specifically, we assign to each 
host halo a collapse epoch Oc, by fitting the mass accretion 
history of the halo, computed as described above in § 13.11 
with the functional form 



Afvir(a) = Mvir(ao) exp(-2ac[ao/a - 1]), 



(3) 



where ao is the scale factor at which we observe the halo 
to have mass Mvir(ao). From the fitted value of Oc, we 
assign each halo a concentration through 



5.125 



(4) 



The mass accretion histories of the subhalos are generally 
more poorly sampled than the accretion histories of hosts 
because subhalos have significantly smaller masses than 
their hosts and are built by correspondingly fewer merg- 
ers above Mmin- It is possible to sample the merger trees 
of halos and smaller subhalos roughly uniformly by low- 
ering the value of Mmin for the smaller subhalos but this 
leads to a rapid increase in the computational expense of 
generating the merger trees and entails calculating merg- 
ers with halos that are far less massive than any halos 
that we expect to host luminous galaxies. Instead, we as- 
sign concentrations to subhalos in the following way. If 
the subhalo mass is greater than lO^Mmin, then we nearly 
always have a well-sampled mass accretion history with 
> 50 merger events and we use the procedure of |W02| as 
outlined above. Oth e rwise , we use the analytic model of 
iBullock et all (|2001bl . EoH hereafter) to compute the mean 
concentration for the mass of each subhalo at the time of 
accretion. We then compute the actual value of Cvir that 
we assign to each subhalo by selecting randomly from a 
log-normal probability distrib ution for cw_with standard 
deviation cr(log(cvir)) = 0.14. iBOll and lWO^ determined 
this to be a good approximation for the statistical scatter 
of Cvir for halos of fixed virial mass in their simulations. 

3.3. Orbital Evolution of Subhalos 

With the accretion histories and the density profiles of 
the host and all accreted subhalos in place, the next step 
is to track the orbital evolution of each accreted system 



to determine how severely each satellite halo is affected by 
dynamical friction and tidal mass loss. 

The first step in tracking the subhalo evolution is to 
assign the initial parameters of the subhalo orbit. Upon 
accretion onto the host, each subhalo is assigned an initial 
orbital energy based on the range of binding energies ob- 
served in numerical simulations. We use the simula tions of 
the MW-like halos Gi, G2, and G3 (IkoI iKGKOl to de- 
termine the initial orbital parameters for subhalos at each 
accretion event. For each accretion event, we located the 
timestep in the simulation at which the most bound parti- 
cle in the subhalo first appeared within the virial radius of 
the main host. We then measured the orbital energy and 
angular momentum of the infalling subhalo. Note that 
because outputs of only 96 epochs were saved for the sim- 
ulation, each infalling subhalo does not have exactly the 
same relation to the host halo at the time that we measure 
these parameters and this may contribute some degree of 
uncertai nty to the ini t ial coii dition s. As in the previou s 
work of iKlvpin et all (|1999a[ ) and IBullock et alt (|2000[ ). 
we find that the initial energy distribution of orbits can 
be well-described by placing each satellite on an initial or- 
bit with an energy equal to the orbital energy of a circular 
orbit of radius i?circ = vRvh , where i?vir is the virial radius 
of the host halo at the time of accretion, and 77 is drawn 
from a uniform distribution on the interval [0.6, 1.0]. 

It is very important to describe the angular momenta of 
satellite orbits correctly in order to estimate the degree to 
which the subhalos are affected by tidal mass loss because 
mass loss is most rapid at pericenter, leading to a very 
strong preference for subhalos on highly e ccentric orbits to 
be disrupted by tides (for examples, see iTavlor &: Babull 
I2OOII l2004aD . We parameterize the initial angular mo- 
menta of subhalos with the orbital circularity e, defined as 
the angular momentum in units of the angular momentum 
of a circular orbit of the same energy, e = J/Jciic- The 
circularity parameter is clearly restricted to the interval 
< e < 1. In Figure [H we show the measured distribution 
of e from the mergers experienced by Gi, G2, and G3 in 
terms of the fraction of in-falling satellites with a given 
initial circularity, d/sat(e)/de. We find the distribution of 
orbital circularities to be well-fit by a one-parameter (3- 
distribution, 



d/,at(e) ^ r(2a) 
de r2(a) ^ 



(5) 



which is valid on the interval < e < 1 and has the simple 
properties of the mean (e) = 1/2, and standard deviation 
cr(e) = l/2\/2a -I- 1. Our best fit of a = 2.22 is shown as 
the solid line in Fig.[T] The distribution of Eq. [S] matches 
the measured mean (e)siin — 0.495 and standard deviation 
o'sim(e) — 0.219 ((Tfit(e) — 0.214) of subhalo orbital circu- 
larities in the simulation quite well. As we were completing 
our work, a s tudy of the init ial orbits of satellite halos was 
pubhshed bv iBensonI (|2004D . the results of which appear 
to be in approximate agreement with our findings. 

The distributions of initial orbits that we have described 
were culled from simulations of three MW-size halos and 
may not be generic. However, over the subhalo mass 
range that we are able to probe, we detect no statistically- 
significant variation in the input distributions with mass. 
Additionally, we observe no significant evolution of these 
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Fig. 1. — Input orbital circularity distribution for infalling sub- 
halos measured from the high-resol ution sim ulations of the Gi, G2, 
and G3 MW-size halos discussed in lKGKO^ . The filled circles show 
the measured distribution for all merging subhalos at the time of 
accretion. The open squares represent the distribution measured 
from only the earliest 25% of all mergers. The solid line shows a 
fit of the distribution of all mergers to the /3-distribution given by 
Equation JSJ with a = 2.22 (see main text for details). 



persion of the host particles by assuming an isotropic 
dispersion tensor; a simple fitti ng form ula for a{r) ob- 
tained in this way is given by IZB03I . We assign the 
Coulomb logarithm ln(A), acco r ding t o the prescription of 
[Hashimoto. Funato. fc Makinol (|2003f ). which they found 
to provide a good match to the results of their iV-body 
experiments for individual orbits: 



In (A) = In 



-M(i?sat/rs), 



(7) 



where r is the radial position of the satellite, i?sat is the 
radius of the satellite (which we identify with the trun- 
cation radius as described below), Ts is the scale radius 
of the NFW profile of the satellite halo, and the function 
liRsat/rs) is an integral of the scattering of background 
particles with impact parameters b < i?sat over the ra- 
dial de i isity p rofile of the subhalo, following the method of 
IWhitd (|1976D . A convenient fitting function for /(i?sat/?'s) 
for NF W profiles is giveii inlZB03l [se e their Eq. (8)]. 

As in iTavlor fc Babull (|200H ) and IZB03I . we estimate 
mass loss due to tidal stripping using a tidal approxima- 
tion. The first step is to estimate the tidal radius rtidc, as 
the radius at which the tidal force from the h ost balance s 
the attractive force of the satellite halo (e.g. lKin"3ll962[ ). 
At each timestep in the subhalo orbit, we estimate the 
mass outside of the tidal radius Afsat(> '"tide), and strip 
mass on a timescale proportional to the inverse of the an- 
gular velocity of the satellite. Specifically, the change 
in mass associated with the satellite at each timestep is 



distributions with redshift. An example of this last point 
is shown in Fig. [Tl where the open squares represent the 
circularity distribution measured for the subhalos with the 
earliest 25% of accretion epochs. The agreement with the 
full distribution is apparent. We proceed by assuming that 
we can set initial conditions by drawing randomly from the 
energy and angular momentum distributions that we have 
described for all merging satellites of all masses at all red- 
shifts. The comparisons with simulations discussed below 
show that this assumption appears to be reasonable. 

We evolve the orbits of satellites by treating them as 
point particles under the influence of the spherical, NFW 
potential of the host halo with an added force term to ap- 
proximate the effect of dynamical friction. As in IZB03I . 
we account for dynamical friction by using a modified 
form of the Chandrasekha r formula (| Chandr asekhadFl 943t 
iBinnev fc TremaindllQSTt ). The approximations of Chan- 
drasekhar result in a frictional force exerted opposite the 
velocity of the satellite of magnitude 



DF 



ln(A) ( erf (a;) - exp(-a;^^ 



where r is the position of the satellite halo relative to the 
center of the host potential, Mgat is the bound mass of 
the satellite, phost{T) is the mass density of the host halo 
at this position, Vsat is the velocity of the satellite, and 
X = Vsat/\/2(T(r)2, where cr{r)'^ is the one-dimensional 
velocity dispersion of particles in the host halo. For the 
NFW host, we estimate the one-dimensional velocity dis- 



<5Afsat = aM,at(> rudc)i2n6t/Lj-^). (8) 

IZB03I motivated the ~ 2Truj^^ choice of timescale by ex- 
amining the typical orbital energ ies of tidally-st ripped 
material in controlled simulations ()Johnstonlll998f ). The 
parameter a is intended to absorb many of the com- 
plicated details of the process of tidal mass loss, which 
is dependent upon the details of subhalo structure (e.g. 
iKazantzidis et al.ll2004af ). and a is the only parameter that 
we fix in order to normalize the model to simulation re- 
sults. We choose a = 3.5 as our normalization and discuss 
this choice further in § 14.11 

Re cent controlled A^-body simulations (jHavashi et al.l 
|2003() and the cosin ological simulations of the Gi, G2, and 
G3 halos (|KGK04| ) indicate that as mass is lost from a 
satellite not only is mass from outside of the tidal radius 
removed, but there is a significant amount of mas s redis- 
trib ution wi thin the tidal radius. iHavashi et al.l l|20Q3[ ) 
and IKGKOJ find that, on average, the maximum circular 
velocity of the stripped subhalos varies with bound mass 
in a remarkably simple way, namely Vmax c>c M^^, with 
7 ~ 1/4 — 1/3 and this scaling has been confirmed in the 
experiments of IKazantzidis et al.l (|2004bl S. Kazantzidis, 
private communication) with much higher resolution. We 

attempt to take approximate account of this redistribution 

1/3 

by enforcing the empirical scaling relation Vinax oc M^^^ . 
For each subhalo, tidal mass loss is most dramatic at 
each pericenter pass. At each subsequent apocenter pass, 
we rescale the density profile of the subhalo by scaling 

1 /3 

Knax oc Mgj^j and set a new scale radius according to the 
relationship between V^ax and Tg for field halos provided 
by the model of iBOll As mass is lost, we assume that the 
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shape of the profile maintains the NFW form [Eq. ([2])] with 
a truncation radius i?sat that we set equal to the radius 
that contains the remaining bound mass, rather than the 
tidal radius, which does not vary monotonically. 

As the final ingredients to our semi-analytic treatment 
of halo substructure, we impose criteria for declaring sub- 
halos to be either completely tidally disrupted or centrally 
merged with the host halo, and therefore no longer identifi- 
able as distinct substructure. We consider a subhalo com- 
pletely tidally disrupted if, after any episode of tidal mass 
loss, the remaining bound mass of the satellite is less than 
the mass contained within its scale radius, Msat(< rg). 
We declare a subhalo centrally merged with its host halo 
if the apocenter of the subhalo orbit becomes smaller than 
f merge — 5 kpc due to dynamical friction. This is a rather 
conservative central merger criterion, and in practice it is 
rarely important in our model because subhalos are always 
severely affected by tides long before their apocentric dis- 
tances fall below rmerge, so that they drop out of any mass 
or 14iax cut that is designed to select subhalos that may 
potentially host luminous galaxies. The purpose of this 
criterion is mostly one of pragmatism, to eliminate the in- 
tegration of very tightly bound orbits which require very 
small time steps. 

3.4. The Subhalos of Subhalos 

When each subhalo merges with a larger host, the sub- 
halo itself may be the host of still smaller satellite halos. 
These subhalos-of-subhalos (SOShalos) can be an impor- 
tant contribution to the total subhalo population, espe- 
cially in the case of a very massive host halo, such as when 
a halo typical of a small group merges into the halo of a 
large cluster. All of the interactions of the SOShalos can- 
not be accounted for in detail without leading to a very 
rapid increase in the computational expense of the calcula- 
tion, mitigating the utility of the semi-analytic approach, 
the raison d'etre of which is speed and simplicity, and so 
so me simplifying a s sumpti ons must be made. 

iTavlor fc Babuil (|2004al) treat SOShalos using an effi- 
cient averaging scheme that is based on the assumption 
of self-similarity of the subhalo population of the primary 
host halo and the SOShalo populations of subhalos. At the 
level of the main host, they determine the mean amount of 
mass lost by subhalos during their evolution and the frac- 
tion of subhalos with a radial position greater than the 
radius containing this amount mass fs- At each merger 
event, they then strip a fraction fs of the SOShalos asso- 
ciated with the incoming subhalo, starting with the most 
recently acquired SOShalos (i.e. last in, first out), and 
place them on correlated orbits within the main host. The 
remaining SOShalos, the fraction 1 — /g of SOShalos that 
were earliest acquired by the incoming subhalo, are then 
considered too tightly bound to be stripped from the sub- 
halo by the main host and are added to the subhalo system. 

We adopt an alternative approximation scheme that is 
more computationally expensive but does not rely upon 
the assumptions of self-similarity, the "last in first out" 
stripping of SOShalos, or the computation of a mean frac- 
tion of stripped SOShalos. However, we make several in- 
dependent simplifying assumptions as follows. We treat 
SOShalos by considering their evolution through three dis- 
tinct phases punctuated by abrupt transitions. First, be- 



fore a subhalo is accreted onto a larger host, it plays the 
role of a host halo and we simply evolve the orbits of any 
SOShalos as described in ? I3.3I Accordingly, SOShalos can 
be significantly dynamically evolved and can even fall be- 
low our minimum mass or Vmax thresholds prior to being 
acquired by the main host halo. Second, upon accretion of 
a subhalo onto a host, we treat the SOShalos of the sub- 
halo, if they are present and still above the minimum halo 
mass of interest, by continuing to integrate their orbits 
within the subhalo potential, neglecting their interactions 
with the main host halo. Third, we consider the possibil- 
ity that a SOShalo may be stripped from the subhalo and 
deposited in the main host. This represents a transition 
where the dynamics of the SOShalo go from being domi- 
nated by the subhalo to being dominated by the main host. 
During the evolution of the subhalo in the main host, if 
the tidal radius of the subhalo becomes smaller than the 
radial position of a SOShalo, we declare the SOShalo to 
be stripped off of the subhalo. We then remove it from 
the subhalo, place it in the potential of the main host, and 
continue integrating its orbit until the final epoch (usually 
z = 0). 

When we place a stripped SOShalo within the potential 
of the main host, we must assign it an orbit that is closely 
associate d with the orbit o f the subhalo from which it was 
stripped. iJohnstoiil (|1998f ) found that the typical specific 
energy of tidally stripped particles is set by the change in 
the host halo potential on the length scale of the orbiting 
satellite, Ai^ts = '^tided$host('')/dr. We set the stripped 
SOShalo on an orbit with specific energy equal to the or- 
bital energy of its parent subhalo in the main host ±Ai?ts, 
with the sign of the shift in orbital energy chosen randomly 
to correspond to the leading and trailing tails of the tidal 
debris. 



4. MODEL RESULTS AND COMPARISONS WITH 
NUMERICAL SIMULATIONS 

In this Section, we present the results of the substruc- 
ture model described in § [3] along with a comprehensive 
comparison with the results of high-resolution numerical 
simulations over the mass range where the two approaches 
can be compared robustly. 

4.1. Subhalo Velocity and Mass Functions 

As in |ZB03| . the particular algorithm described in § [3] 
made use of inputs from numerical simulations and, in 
particular, was designed to approximately reproduce the 
cumulative velocity fun ction (CV F) A'sat(> ^ax), of sim- 
ulated MW-size halos (|KGK04D . The CVF is defined as 
the number of subhalos above a given threshold in max- 
imum circular velocity Knax, as a function of T^nax- Fig- 
ure [2] shows model predictions for CVFs compared to the 
results of numerical simulations. In all cases, we refer to 
subhalo maximum circular velocities V^ax; that are scaled 
by units of the maximum circular velocity of the host halo 
Knax* ; Order to approximately scale out the dependence 
of the velocity function on the size of the host halo and 
facilitate the comparison of velocity functions of hosts of 
various sizes. The solid line in Fig.[2]shows the mean CVF 
from 100 model realizations of a host halo with virial mass 
Mhost = 10^^-^ /i-^Mq in a standard ACDM cosmology 
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with ^ I - ^ 0.3, h ^ 0.7, and as = 0.9. The 
error bars on this result represent the 1 — cr realization-to- 
realization scatter about the mean value measured from 
these 100 model realizations. The dotted lines show the 
measured cumulative velocity functions of the three sim- 
ulated MW-size halos Gi , G2, and G3 introduced in § O 
The dotted lines are truncated at the value of T/^ax where 
the simulation halo catalogs become incomplete due to fi- 
nite resolution. The good agreement between the model 
CVF and the CVFs of the simulated MW-size halos is 
apparent and was achieved by setting the mass l oss pa- 
rameter a = 3.5 in Eq. [51 We emphasize that, as in [ZB03l . 
our semi-analytic model was normalized to this statistic 
only and as such, it is now possible to make predictions 
for other statistics, and to diagnose and assess the general 
applicability of the model by comparing these predictions 
to the results of numerical simulations. 

Figure [3 also shows the first prediction of the semi- 
analytic model compared to the result of an A^-body 
simulation. The dot-dashed line shows the CVF of a 
A^host = lO^**'^ h~^MQ cluster-size host halo in the same 
"concordance" cosmology. For clarity, the scatter among 
these realizations is not shown, but it is similar to the scat- 
ter on the CVFs of the MW-size halos. The dashed line 
is the CVF of a simulated cluster-size halo of the same 
mass. The predictions of the model and the simulation 
are in good agreement as the simulated halo is well within 
the 1 — a scatter about the mean CVF of the model. No- 
tice also that the model predicts that the CVF should not 
scale in a precisely self-similar manner with the size of the 
host. The model results in a CVF that is slightly higher 
at higher values of ■ We return to this point shortly. 

As a further test of our model, we now turn to the cu- 
mulative mass function (CMF) A'sat(> -A^sat) of subhalos. 
As with the CVF, we refer to subhalo masses Mgat in units 
of the host halo mass Mhost in order to scale out the gross 
dependence of iVsat(> Afsat) upon the mass of the host 
halo. We show CMFs in Figure [3] Again, the agreement 
between the results of the model and the numerical sim- 
ulations is good, lending confidence in the ability of the 
model to account approximately for the dominant phys- 
ical effects that determine subhalo properties, and to be 
applied more generally. Note also that, just as we saw 
for the CVF in Fig. H the CMFs in Fig. [3] show that 
the number of subhalos of a given size does not scale self- 
similarly with the size of the host. Larger host halos tend 
to have more substructure at a fixed value of Mgat /Mhost 
(or VJ^^X/ V^°l^)^ in agreeme nt with the recent numerical 
results of Gao et al. ( 2004a') and the simplif ied analytic 
considerations of .van den Bosch et al.l (|2004bO . 

4.2. The Radial Distribution of Subhalos 

The radial distribution of subhalos within their hosts 
has received a significant amount of attention re- 
cently and is a necessary ingred i ent o f models of 
galax y clustering statistic s ( SeliakI I2OOOI : ISheth et"al] 
I2OOII: iBerlind fc Weinberg! |2002l ) and in the model- 
ing of the flux-rat io anomalies i n strong lens sys- 



tems (e.g., iDalal fc Kochanek 20"0ll : iBradac et all l2002t 



IChen. Kravtsov. fc Keetoiil 12003^ At least three groups 
have recently performed convergence studies that demon- 
strate that the radial distribution of subhalos of all 
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Fig. 2. — The mean cumulative velocity functions A^sat (> Vmax), 
of satellite halos. The solid line depicts the mean value of the 
number of subhalos with a maximum circular velocity greater than 
y^axi ^ ^ function of VJ^^x predicted by 100 model realizations 
of a Mhoat = 10^^'^ h~ Mq host halo. Velocities are scaled in 
units of the maximum circular velocity of the host halo V^°J* . The 
errorbars depict the I — a scatter among the 100 model realiza- 
tions. The dotted lines show the cumulative velocity functions of 
the three simulated MW-size halos. The dashed line shows the cu- 
mulative velocity fun ction of t he simulated cluster-size halo (CL2 
HR in Tasitsiomi et al.l 120041 ). The dot-dashed line is the mean 
predicted cumulative velocity function from 100 model realizations 
of a Mhost = lO""^*'^ h~^yiQ halo. In the interest of clarity, the 
realization-to-realization scatter is not shown for this case but is 
similar to the scatter for the MW-size model realizations. 



sizes is significantly less central ly-concentrated than 
the distribution of dark matter (jPiernand et al.l |2004 
iGao et ani2004aHNagai fc Kravts ov'200 5D, confirming sev- 



eral e arlier results (iGhigna et al.lil998t. 12000 ; IColm et all 
19991: iSpringel et all |2001D . iDiemand et al i ([200^ and 



Gao et al.l (l2n04bh stressed that the observed spatial distri- 
butio n of galaxies in clusters ((Carl berg. Yec, & EUingsoni 
1997t Ivan der Marel et all l2000t lETn. Mohr. fc' Stanford 
2004[) appears to be more centrally-concentrated and 
not as strongly anti-biased with respect to the dark 
matter as the distribution of subhalo s in dissipation - 
less iV-body simulat i ons (se e also Springel et al.l 1200 1[ ). 
Ivan den Bosch et"aLl (|2004a[ ) analyzed satellite galax- 
ies in the Two-De gree Field Galaxy Redshift Survey 
(|Colless et al.l I2OOII ) finding preliminary indications that 
observed galaxies do follow a centrally-concentrated, 
NFW-like radial distribution over a wide range of host halo 
masses, though they cannot rule out shallow er distribu- 
tions because of incompleteness of close pairs (jCole et al.l 
I2001D and more detailed studies are required to confirm 

this. 

iTavlor Silk, fc Babull (|200l and iTavlor fc Babul 
(j2004bl lcD used semi-analytic methods to argue that 
the distribution of subhalos should be more centrally- 
concentrated than that observed in simulations and 
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Fig. 3. — The mean cumulative mass function (CMF) of sub- 
halos as a function of subhalo mass relative to the mass of the 
host halo, Msat /JV/jjost . The solid line shows the mean value of 
{Ns^,t) (> Msat/Mhost) for a Mhost = lO^^'^ H-'^Mq host computed 
from 100 realizations of the dynamical model of § [3] The errorbars 
depict the I — a scatter among the 100 realizations. The dotted 
lines show t he C MFs of the t hree simulated halos of similar mass 
described in lKOll and lKGK04l . The dot- dashed line shows the mean 
CMF of a Mhost = 10"'2 H-'^Mq halo predicted by 100 model 
realizations (for clarity, the scatter is not shown in this case but is 
similar to that of the smaller host). The dashed line shows the CMF 
measured from a high-resolution simul ation of a c luster-size halo of 
a similar mass (halo CL2 HR in iTasitsiomi et al.ll2004il . 



simulations may still be subject to the problem 
overmerg inK" in the central regions of halos (see 
The r es ults of the recent stud- 
2004f ). iGao et all (l2004aD . and 



that 
of " 

iKlvpin et al 



T999bh. 



ies by [ Diemand et al 



iNagai &: Kravtso v (200 ^ cast doubt upon the proposal 
of lTavlor et'aLrft2003. ). .Nagai fc Kravtsov! (|2005D suggest 
that the radial distribution of galaxies selected on lumi- 
nosity cannot be compared directly with that of subhalos 
selected by mass because galaxies do not likely lose 
luminosity at the same rate that they lose dark matter 
while they orbit inside a larger host halo. Nevertheless, 
it is important to pursue further tests of this issue. 
The growing importance of the radial distributions of 
subhalos and various subsets of subhalos with specific 
properties, coupled with these recent theoretical studies, 
provide interesting motivations for comparing the spatial 
distributions of subhalos in our model with that observed 
in numerical simulations. 

In Figure [H we compare the results of our model pre- 
dictions for the radial distributions of subhalos with the 
results of the Gi, G2, and G3 MW-size halo simulations 
and the high-re s olutio n cluster simulation discussed in 
ITasitsiomi et all (|2004[ ). In Figure H we plot the cumu- 
lative number of subhalos at a distance less than R away 
from the center of the host halo Nsa,t{< i?), as a function 
of R and normalize by the total number of subhalos within 
the virial radius of the halo iVtotai (Fig.[2]and Fig.[3]already 



show that the values of iVtotai are in agreement). In the left 
hand panel, we compare the radial distribution of subhalos 
within MW-size host halos. We choose the selection crite- 
rion V^tx/^nax ^ 0-07 in order to restrict our results to 
subhalos that are sufhciently large as to be well-resolved in 
the simulation and fairly insensitive to numerical effects. 
The dot-dashed line shows the distribution of dark mat- 
ter for an NFW profile with a concentration parameter of 
Cvir = 12.6. This is a t ypica l value of Cvir for a MW-size 
halo in this cosmology ( BOll ) and is the mean value of the 
concentration of the 100 hosts in the semi-analytic model. 
The three simulated MW-size halos Gi, G2, and G3 have 
best-fit concentration parameters of Cvir — 9.5, 12.5, and 
14.5 respectively. Similarly, in the right hand panel of 
Fig. [U we compare the radial distribution of subhalos of 
the simulated cluster halo with the distribution of subhalos 
in 100 model realizations of a Mhost = lO"'^'*'^ h~^MQ host 
halo. The maximum circular velocity threshold in this case 
is t^max/^max* > 0-05. Oncc again, the dot-dashed hne rep- 
resents the distribution of dark matter for an NFW profile 
with Cvir = 6.1, which is typical for a halo of this mass 
and is close to the mean concentration of the 100 hosts in 
the semi-analytic model. The simulated cluster halo has 
a best- fitting concentration of Cvir = 8.5, which is roughly 
^ 1 — (T higher than the mean concentration at this mass 
according to the halo-to-halo scatter observed in simula- 
tions (BOl). For the simulated subhalo distributions, we 
computed profiles by stacking the last five simulation out- 
puts (a = 1.00,0.99,0.98,0.97,0.96 for Gi, G2, and G3 
and a = 1.000,0.995,0.990,0.985,0.980 for the cluster) in 
order to overcome the noise in the measurement. 

Several interesting features are evident in Figure [31 
First, the anti-bias of the subhalo population with respect 
to the dark matter is evident in both panels of Fig. 2] for 
both the simulated halos and the model halos. A MW-size 
halo typically has ^ 35% of its mass within ~ 0.2i?vir, yet 
it has only ~ 5% of its satellite halos located this close to 
its center. Second, notice the excellent agreement between 
the simulation and model subhalo profiles in both panels. 
Again, this agreement is a non-trivial confirmation of the 
success and more general applicability of the semi-analytic 
model. After designing the model to reproduce the num- 
ber of subhalos in a MW-size host, the model naturally 
recovers the correct radial distribution of satellites within 
both a MW-size host halo and a more massive cluster-size 
halo. 

The agreement between the model subhalo distribution 
and the distribution of subhalos in the simulations is an 
encouraging sign that the semi-analytic model captures 
the dominant physical processes that determine the de- 
mographics of CDM subhalos. Moreover, this indicates 
that the model of § [3] results in a spatial di stribution of 
subha los t hat is in conflict with th e model of lTavlor et al.l 
([2003') and lTavlor fc Babull (|2004b[ ) , who find a much more 
centrally-concentrated distribution of subhalos. The radial 
distributions of subhalos in our model support the conclu- 
sion that the anti-bias of subhalos with respect to dark 
matter is a physical effect rather than the result of numer- 
ical overmerging, yet it will likely require further study 
to provide a definitive answer to this question. However, 
our results indicate that the disagreement between semi- 
analytic models and simulation results may depend upon 
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0.1 1 0.1 1 

Fig. 4. — Left Panel: The cumulative radial distribution of satellite halos in MW-size host halos. We plot the fraction of subhalos within 
a given distance from the host halo center Afsat(< ^i)/Aftotal a- function of distance R, in units of the halo virial radius, R/R^iT- The 
solid line shows the mean cumulative radial distribution of subhalos with V^^x ^ 0.07Vj5J°J for Mhost = lO^'^'^ h~^MQ computed from 100 
realizations of the semi-analytic substructure model described in § [3] The error bars show the la scatter in the 100 model realizations at 
only four radii for clarity. The dotted lines show the cumulative radial subhalo distributions for the Gi, G2, and G3 simulated MW-size 
halos. The da shed line represents the distribution of dark matter for an NFW halo with concentration parameter Cvir = 12.6, typical of a 
MW-size halo lIBOlf) . The simulated halos have best-fitting NFW concentrations of c^ir ~ 9.5, 12.5, and 14.5. Right Panel: The cumulative 
radial distribution of satellite halos in cluster-size host halos. The line types are the same as in the right panel. In this case, the host mass 
is Mhost = 10^*'^ h~^MQ and all subhalos with V^^x ^ 0-05V^°J* are plotted. The dashed line shows the distribution of dark matter for 
an NFW profile with concentration Cvir = 6.1, which is typical for this host mass. The best-fitting concentration for the simulated halo is 
Cvir — 8.5. In both panels, the simulated radial distributions were computed by stacking the last five simulation outputs for each object. 



the specifics of the model implementation and is not a 
generic feature of such models. 

4.3. The Halo Occupation Distribution of Subhalos 

We now turn to a more comprehensive comparison of 
our semi-analytic model to the predictions of cosmological 
A^-body simulations. To do this, we compare the results 
of our model to the subhalo populations in the ACDMgo 
simulation studied by KraytsOY et al. ( 200^ . Recall that 
the ACDMgo simulation was performed for a ACDM cos- 
mological model with as = 0.75. This power spectrum 
normalization is lower than the normalization of erg = 0.9 
used in the simulations of the Gi , G2 , and G3 halos as well 
as the CL2 cluster halo to which we compared our model 
in § 14.11 and § 14.21 We therefore run our model again us- 
ing cTg = 0.75. The parts of the model that are directly 
affected by this change are the EPS merger trees of halos 
and their mass density profiles. Comparing the ACDMgo 
simulation results with the model at this different power 
spectrum normalization and at different halo masses than 
the CTg = 0.9, Mhost = 10^^-^ h~^MQ simulations that we 
used to tune the model is not a trivial exercise and tests 
the general applicability of the model. 

Our first comparison is shown in Figure O In this plot, 
we show the mean number of subhalos, (iVgat), as a func- 
tion of host halo mass, for four different T4iax thresholds: 
Vmax > 100, 150, 200, and 250 kms^^. The lines show 



the result of 1000 model realizations in each host mass bin 
from log(Mhost/ /i-^Mq) = 11.0 to log(Mhost/ H-^Mq) = 
15.0 with bins evenly-spaced by A (log Mhost) = 0.1. The 
squares and triangles represent the (A^sat) measured from 
the halos in the ACDMgo simulation. 

Figure [5] shows a rather remarkable result: our model of 
subhalo populations matches the results of the simulations 
very well at each host mass and each Vmax threshold. At 
a single, fixed value of Mhost, the four model or simula- 
tion points associated with that value correspond to four 
points on the CVF for host halos of that mass. As such, 
Figure [5] demonstrates that, after culling the input dis- 
tributions for impact energies and impact parameters of 
subhalos at accretion from a cosmological A'^-body simu- 
lation of the formation of three MW-size halos, and fixing 
our one-parameter model to match the amplitude of the 
CVF of subhalos of a MW-size host, the model then suc- 
cessfully reproduces the CVFs of subhalos over a range 
of roughly three orders of magnitude in host mass in a 
cosmology with a different power spectrum normalization. 

As a further diagnostic, we can go beyond the first mo- 
ment, or mean, of the occupation distribution of subhalos 
and explore the shape of the probability distribution for 
hosting iVsat satellites at a given host mass P(iVsat|Mhost), 
by comparing the second moment of the distribution. The 
higher order moments of the subhalo distribution are im- 
portant for studies of the iV-point statistics of subhalo (or 
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Fig. 5. — The mean subhalo occupation number as a func- 
tion of host halo mass, for several peak circular velocity thresh- 
olds: Vmax > 100, 150, 200, and 250kms-i. Squares and tri- 
angles show results of the ACDMgo A^-body simulation, while the 
results of our semi-analytic model described in § |3] is shown by the 
solid curves. The model results are based on 1000 Monte Carlo 
realizations in each host mass bin, with host mass bins spaced by 
A(log Mjjogt) = 0.1. For both the simulation and the model results, 
the error bars represent the estimated error on the mean value of 
A^sat in each mass bin. 




Fig. 6. — The square root of the second moment of the halo occu- 
pation distribution of subhalos (Afsat(Afsat — l))^^^i ^ a function of 
A^host ■ The solid lines are the mean of 1000 realizations of the semi- 
analytic model of § [3] The squares and triangles show the results 
for the subhalo populations measured in the ACDMgo cosmological 
Af-body simulation. From top to bottom the three sets of lines and 
points are for three different subhalo Vmax thresholds; Vmax > 100, 
150, and 200kms~^ respectively. 



galaxy) populations where N >2. In Figure [51 we show 
results for the second moment of the subhalo occupation 
distribution (iVsat(iVsat — 1)), as a function of host mass. 
As in Fig. [51 the lines represent the model results, the 
points represent the simulation results and the error bars 
represent the error on (iVsat(-^sat — !))■ We show sub- 
halo populations above three cuts in Knax, subhalos with 
Knax > 100, 150, 200 kms~^. The agreement between the 
model and the simulations is good for the second moment 
of the distribution of subhalo populations for the highest 
two Finax cuts. Howcvcr, a notable feature of Fig. [51 is the 
discrepancy between the model results and the simulation 
results at low values of Mhost or, equivalently low values 
of {Nsat{Nsat — 1)), in the lowest Knax sample, where the 
model overpredicts the second moment of the distribution 
b y as much as ^ 4 0%. 

iKravtsov et al.l (|2004f ) found the second moment of the 
occupation distribution of subhalos to be consistent with 
that of a Poisson distribution in the host mass regime 
where (A^sat) ^ 0.1. In this case, the probability distri- 
bution function for a host halo of mass Mhost to con- 
tain a particular number of satellite halos is completely 
specified by the mean of the distribution because the 
higher order moments are related to the mean through 
(A^(7V- 1) ...{N-i)) = (7V)*+^ In Figure[71 we compare 
the occupation distribution of subhalos predicted by our 
semi-analytic substructure model to a Poisson distribution 
directly by plotting the quantity 



a EE (iVsat(iVsat - 1)) V (^^.at) • (9) 

1 /2 

For a Poisson distribution, (Afsat(A^sat - 1)) = (A^gat) 
so that a = 1, while for distributions that are narrower 
(sub-Poisson) or broader (super-Poisson) than a Poisson 
distribution, a < 1 and a > 1 respectively. Figure [3 
shows a as a function of Mhost for subhalos in the semi- 
analytic model and the ACDMgo simulation for subhalos 
with Vmax > 100 kms^^, and Vmax > 200 kms^^ in the left 
and right panels respectively. Figure [3 shows that for the 
Knax > 100 kms~^ threshold, the semi-analytic model pro- 
duces a probability distribution for the number of satellites 
at a given host mass that is significantly super-Poisson in 
the host mass range 10^^^ Mhost/ H'^Mq^ lO^'^, corre- 
sponding to a mean number of satellites 0.1^ (A"sat) ^ 1- 
In this range, the model over-predicts a for the occupation 
distribution of subhalos by as much as ^ 25% with respect 
to the simulated occupation distribution of subhalos. The 
simulated occupation distribution of subhalos is consistent 
with the Poisson value of a = 1 for A/hostJc] 10^^ /i^^Mq 
((A^sat) ^ 0.1). At higher Kiax thresholds, the magnitude 
of the discrepancy decreases, as can be seen in the right 
panel of Figure[71 but qualitatively, we find that the super- 
Poisson "bump" at (Nsat) 0.1 — 1 persists at all T4iax 
thresholds. The over- prediction of the second moment of 
the distribution relative to the simulation result appears 
to be a fairly general feature of our semi- analytic model. 
The model of lvan den Bosch et al.l ()2004b[ ) yields a similar 
overestimate of the second moment despite the fact that 
their mass loss model is quite different from the model 
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Fig. 7. — The width of the halo occupation distribution of subhalos compared to the width of the Poisson distribution. We plot the 
quantity a = (A'sat(Afsat — 1))"'^''^ / (A^sat) as a function of host mass for two circular velocity thresholds: Vmax > lOOkms"^ (/e/t panel) 
and Vmax > 200kms~^ {right panel). The solid lines correspond to the results of the analytic model of §|3]and the squares correspond to 
simulation results. The horizontal dotted line marks the value a = 1 which corresponds to the Poisson distribution. 



that we present, and is based on a simplified averaging of 
subhalo mass loss (F. C. van den Bosch, private commu- 
nication). Coupling this with the fact that EPS formation 
times seem to exhibit a larger scatter than the f ormation 
times of simulated halos, as pointed out bv I W02l leads us 
to speculate that the origin of this discrepancy is a fun- 
damental shortcoming of the standard EPS formalism for 
generating halo merger trees. Further investigation is nec- 
essary to test this hypothesis. 

As a further test of the model, we generalize this com- 
parison and examine the subhalo occupancy predicted by 
the model as a function of redshift. First, Figure [5] shows 
the mean number of subhalos with V^nax > 100 km s^^ as 
a function of host mass at redshifts z = 0, 1, and 3. We 
illustrate the redshift dependence using only this one Knax 
threshold at high redshift because the number of large ha- 
los decreases rapidly with redshift, causing the statistics 
to be rather poor at larger thresholds. At z = the model 
points represent the same 1000 realizations per mass bin 
as shown in Fig. [5] At redshifts z = 1 and z = 3, the 
model results are based on 100 realizations per mass bin. 
Again, the agreement between the model and the ACDMgo 
simulation is rather impressive. The simple model of halo 
substructure that we developed in §[31 predicts a mean sub- 
halo occupation number at high redshift that is in good 
agreement with the predictions of the full numerical sim- 
ulation. This result also demonstrates the clear trend for 
CDM halos to host a larger amount of substructure at 
higher redshifts. Here we select halos at a fixed Vmax cut; 
using a sample selected at a fixed number density, yields 
the same trend, though slightly less pronounced. 

The high-redshift agreement is not limited to the mean 
occupation. Figure [9] shows the corresponding second mo- 
ment of the occupation distribution of subhalos as a func- 
tion of redshift. The results are for the same model realiza- 
tions used to generate Fig. [51 The figure shows that unlike 
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Fig. 8. — The redshift dependence of the mean of subhalo occu- 
pation. The lines show the predictions of the semi-analytic model 
at redshifts z = (bottom), 2 = 1 (middle), and z = 3 (top). 
The model results are based on 1000 Monte Carlo realizations at 
each host mass bin at z = 0, and 100 realizations at redshifts 
z = 1,3. The filled squares and open triangles rep resent the same 
quantity measured from the simulations described in lKravtsov et ah! 
(2004). For both the model predictions and the predictions of the 
simulations, the errorbars represent the estimated error on (iVsat). 
Note that we plot the subhalo occupation at only one threshold 
(Vmax > lOOkms"^) shown in Fig. [5] because large halos become 
increasingly rare with increasing redshift causing the statistics to be 
poor at high redshifts and large thresholds. 
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Fig. 9. — The redshift dependence of the second moment of the 
subhalo occupation distribution, (Nan.tiNaa.t — ^))^^'^ ■ The lines show 
the results of the semi-analytic substructure model for subhalos with 
Vmax > 100 km s"-*^ at redshifts 2 = {bottom), 2 = 1 (middle), and 
2 = 3 (top). The model results are based on 1000 model realizations 
at each host mass bin at 2 = 0, and 100 realizations at redshifts 
2 = 1,3. The filled squares and open triangles represent the same 
quantity measured from the ACDMgo simulation. 



at z = 0, the high redshift second moment model results 
are in good agreement with the ACDMgo simulation over 
a wide range of host halo masses. As redshift increases at 
a fixed host mass bin, the masses of the host halos in that 
bin become increasingly large with respect to the typical 
mass of collapsing objects Mi,. So as redshift increases 
at fixed mass, we are examining increasingly rare objects. 
The results in Fig. [9] illustrate the same qualitative trend 
that is seen in Fig. [SI namely, that the model result for the 
second moment appears to be more like the simulation re- 
sults as we examine increasingly rare host halos. 

4.4. Accretion vs. Destruction 

Although the first detailed studies of halo substructure 
showed that subhalo populations may scale in a simple, 
near ly self-similar way with the size of the the host halo 
(e.g. iMoore et al.lll9 9ff) there is now evidence to the con- 
trary. Recently, [Gao et al.i (|2004a|) reported the mea- 
surement of a deviation from the self-si milar scaling of 
subhalo populations in simulated halos (jPiemand et al.l 
|2004 hinted at this but did not have a large enough 
host halo sample to make a significant detection) while 
Ivan den Bosch et alj (|2004bl ) found a similar trend in sub- 
halo abundance using analytic arguments. To be spe- 
cific, iGao et al] ([20043) found that the number of sub- 
halos with mass scaled to a fixed fraction of the host 
halo mass d-/Vsat/d(Msat/Mhost), increases with host mass 
as - MO^i^ in the regime M^at/Mhost < 1- In Fig- 
ure [21 we show that the results of our semi-analytic model 
predict a very similar trend with host halo mass and, 
as we have a large number of realizations of host halos 



in each mass bin, the trend that we compute is signifi- 
cant. O ur model is in good agreement with the lGao et al.l 
(|2004aD result, yielding a subhalo abundance that scales 
as a power-law in the scaled satellite mass (Msat/Mhost), 
with a weakly host mass-dependent normalization. We 
find that dA^sat/d(M,at/Mhost) oc Mh^ost(^^>5at/A4ost)"'' in 
the regime Mgat/Mhost < 1, with v « 0.08 and ^ « 1.88. 

The semi-analytic model also provides a simple frame- 
work for understanding the nature of this trend. The 
abundance of subhalos is determined by the constant com- 
petition between subhalo accretion and destruction. In 
Figure [TUl we illustrate the competition between subhalo 
accretion and destruction rates as a function of Afhost for 
different final (z — 0) masses of the host halo. In the stan- 
dard, ACDM cosmological model, the typical amplitude of 
mass density fluctuations is a decreasing function of length 
scale (or an associated mass scale M — AttpmR^/'^ as in 
§ 13. ip . so structure forms hierarchically as small objects 
collapse early and merge to form larger objects. Large ha- 
los observed at redshift z — have thus typically assem- 
bled most of their mass relatively more recently than their 
less m assive counterparts (e.g., see discussions in LC93 and 
IWO2D . Correspondingly, one important manifestation of 
the more recent mass assembly of relatively massive halos 
is that large host halos have typically acquired their satel- 
lites more recently than less massive hosts. This shift in 
the relative accretion time of subhalos is evident in Fig. [TUl 
The rate of satellite halo accretion (shown by the thick, up- 
per lines) for the -Mhost — 10^^'^ /i~^Mq host is strongly 
peaked at ^lookback ~ 10 — 11 Gyr in the past and drops off 
rapidly at more recent times, while for the largest. Coma- 
size host, with Mhost = lO^**'^ h~^MQ host, the accretion 
rate peaks more recently at iiookback ~ 8 — 9 Gyr ago and 
remains nearly constant from this time until the present. 

The thin, lower lines in Fig. [TUl show the accretion rate 
for subhalos that remain above the cut VJL^X > 0.21/^1°^' 
at redshift z = 0. We refer to these subha l os as sur- 
viving subhal os. In agreement with ZB031 [Gao et al] 



([2004aD . and [van den Bosch et al.[ ([2004bD . we find that 
only ^ 4—12% of surviving subhalos are accreted prior to a 
redshift of z = 1 (tiookback ^ 7.6 Gyr). At iiookback ~ 4 Gyr 
(z ~ 0.35), the probability of a subhalo surviving and re- 
maining in the z = sample is roughly 50% and is only 
very weakly dependent upon the mass of the host halo. 
This timescale for removing subhalos from the sample due 
to mass loss and merging is not surprising because the typ- 
ical dynamical times in halos at low redshift are of order 
^3 — 5 Gyr and lends suppo rt to the simpler approach of 
[van den Bosch et al.[ (l2004b) for modeling subhalo popu- 
lations in applications where the level of detail we present 
here is not needed. 

The reason for the deviation from a self-similar scal- 
ing of subhalo populations is evident in Fig. [TUl Subhalos 
in larger hosts have typically been accreted more recently 
and have had less time to be tidally stripped and less time 
for their orbits to decay via dynamical friction and fall to 
the host halo center. The later formation times of more 
massive host halos cause the subhalo populations of host 
halos to deviate from self-similarity simply because the 
balance between typical accretion times and typical de- 
struction times shifts m onotonically in favor of accre tion 
with increasing Mhost- [van den Bosch et al.[ ([2004bD de- 
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veloped a similar interpretation of these results using an 
analytic model of subhalo populations. 

Given the discussion of the trend in the normalization 
of subhalo mass and velocity functions above, it seems 
natural that the trend of subhalo abundance with red- 
shift shown in Figure [5] should have a similar explana- 
tion. Along the same lines, it is also natural to inquire 
about subhalo survival probabilities for satellites accreted 
at various times. We address both of these issues in Fig- 
ure [TT] In the left hand panels of Fig. [TTl we show the 
accretion rates of all subhalos (thick line) and of surviving 
halos (thin line) for our model realizations of host halos of 
mass Mhost = lO^'^'^ h~^MQ at three different redshifts, 
z = 0, 1, and 3. In this figure, we define our subhalo sam- 
ples by a fixed threshold of T4iax > 80kms"^. We refer 
to subhalos that remain above this threshold at the final 
redshift as surviving subhalos. We show accretion rates for 
halos at disparate redshifts for which the typical timescales 
of accretion and the typical dynamical timescales of ha- 
los are quite different from each other. In order to make 
the results for the halos at different redshifts more nearly 
commensurable, in Fig. [11] we assign subhalos an accretion 
time in units of the typical dynamical time of the host halo 
at the time of accretion tiookback/^dyn- The scaling of tdyn 
with ^lookback is responsible for the accretion rates drop- 
ping monotonically with time, as opposed to rising, as they 
do in Fig. [TOl In each panel, the vertical dotted line shows 
the median value of iiookback/^dyn for all accreted subha- 
los. The results are computed from 1000 realizations of 
the model of § [3l 

Notice first in Fig. [Til the relative shift in the lookback 
time to accretion as a function of the redshift at which 
we observe the host halo. For the z — host halos, the 
median subhalo accretion time is roughly four halo dy- 
namical times, while for the z = 3 hosts, the median halo 
accretion time is only about one-half of a halo dynamical 
time. The reason for this shift is essentially the same as 
the reason for the shift in accretion times as a function 
of host halo mass: large halos are increasingly rarer ob- 
jects at high redshift than they are at low redshift and are 
characterized by more recent mass assembly. 

In the right-hand panels of Fig. [Til we explicitly 
show the fraction of surviving subhalos /surv accreted at 
^lookback, as a function of ^lookback in units of the typical dy- 
namical time of the host. In the upper right of each panel, 
we also show the total, integrated fraction of all accreted 
subhalos that survive until the final redshift, Fgurv- In all 
cases, subhalos very rarely survive for more than '--^ 3 — 4 
dynamical times after they are accreted onto their host 
halos. The typical timescale for destruction is roughly a 
dynamical time, as approximately 50% of all satellite ha- 
los fall below the sample Vmax threshold after roughly one 
dynamical time in each case. However, notice that in the 
hosts observed at a final redshift of z = 0, the typical 
time for a subhalo to be accreted is ^ 3.8idyn, while for 
the host observed at z — 3, the typical subhalo accre- 
tion time is only ~ 0.6tdyn- This shift in the accretion 
times relative to the timescale that is relevant for remov- 
ing a halo from the threshold sample is dramatic and is 
responsible for the increase in the overall subhalo survival 
fraction as i^surv — 0.16, 0.22, and 0.26 for the redshifts 
z — 0, 1, and 3 halos respectively, and it is the reason 



for the systematic increase in the abundance of subhalos 
in fixed mass hosts with increasing redshift illustrated in 

Fig.m 

4.5. Host Properties and Subhalo Populations 

The discussion of the previous section has a natural ex- 
tension. As can be seen in Fig. [21 there is a sizable scatter 
in the number of subhalos of a particular size in a host of 
fixed mass. In the preceding paragraphs, we demonstrated 
that more massive host halos typically contain more sub- 
structure simply because they assembled their mass more 
recently and their subunits have had less time to evolve dy- 
namically. A natural implication of this is that the scatter 
in the amount of halo substructure at fixed host halo mass 
may be largely determined by the variety of mass accre- 
tion histories of halos at fixed mass: halos that acquired 
their mass more recently should have more substructure. 

A convenient quantity that can be used to describe the 
variety of mass assembl y hist ories of halos is the collapse 
epoch Qc, introduced bv lW02l and defined in Eq. ^ above. 
As we described in 13.21 for each realization of our model, 
we determine the collapse epoch of the host halo by fitting 
Eq. ([3|) to its mass assembly history. 

In the left panel of Figure [T^ we show the correlation 
between Oc and satellite number predicted by our model 
for five different values of host mass. For each Mhost value, 
we plot the number of satellites with V^^^ > 0.15V^°^, 
in units of the mean number of satellites for hosts of 
this mass, as a function of Oc in units of the mean col- 
lapse epoch (oc), for hosts of this mass. The thick, 
solid, central line corresponds to the mean relation for 
Mhost = 10^^-^ H'^Mq hosts and the upper and lower 
solid lines show the scatter in this relation measured from 
1000 model realizations. The remaining lines show the 
mean relation between collapse epoch and satellite num- 
ber for the four different host halo masses indicated in the 
Figure caption. The strong correlation between collapse 
epoch and satellite number is apparent in Fig. 1121 with 
host halos that have an average collapse epoch hosting 
an average number of satellite halos. We see no evidence 
that this relation changes significantly as a function of host 
mass. Additionally, at fixed Mhost and Oc, there remains 
a significant amount of scatter in the number of satellite 
halos, indicating that there are other important physical 
ingredients that affect the amount of halo substructure. 

In the context of our model, two possible sources of this 
scatter are: (1) the distribution of orbital parameters and 
subhalo concentrations for a fixed host halo accretion his- 
tory; (2) the variety of accretion histories that result in the 
same best-fit Uc- By constructing 100 model realizations 
using a single host mass accretion history in several mass 
bins, we estimate that approximately half of the scatter 
can be attributed to (1) above. 

As we mentioned in § 13.21 the concentrations of halos 
correlate strongly with their mass accretion histories. In 
fact, we used this correlation to set the concentrati ons of 
our host halos according to Eq. (|3|), as prescribed bv lW02l 
The implication is that the fundamental correlation be- 
tween mass accretion history and satellite abundance im- 
plies a correlation between host Cvir and A'sat, with A'sat a 
decreasing function of host halo concentration. In the con- 
text of our model, this correlation must be present given 
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Fig. 10. — The distributions of subhalo accretion and destruction times as a function of the final 2 = host halo mass Mjjost- In each 
panel, we show accretion rates for the different final host masses listed in the top left of the panel. The thick solid upper lines indicate the 
number of all subhalos that are accreted with an initial maximum circular velocity > 0.2V^°^ as a function of the lookback time to the 

time of accretion (i.e. time since accretion) in bins of width Atiooi^back = 0.5 Gyr. The thin solid lower lines show the number of subhalos 
accreted in each lookback time bin, which remain above the same threshold of V^^x at z = 0, after the various destructive effects have been 
accounted for. The number of subhalos in each bin is shown in units of the total number of accreted subhalos Ntot- 



the correlation between ac and A^sat shown in the left panel 
of Fig. [T2] and the fact that we assign concentrations via 
Eq. ([4]). We illustrate the model results for the Cvir-./Vsat 
relation in the right hand panel of Fig. [121 Again, the 
different lines correspond to the model predictions at vari- 
ous host halo masses. The strong correlation between halo 
concentration and satellite number is clear. For this sub- 
halo selection criterion a power-law fit to the semi-analytic 
model result yields iV^at oc c~-'^ with a « 0.86 ± 0.07. 

It is interesting to test the extent to which the simulated 
halos exhibit this correlation between concentration and 
subhalo number. The stars in Fig. ll2l show the relationship 
between TVgat and Cvir for host halos in the ACDMgo simu- 
lation. The values for the simulation halos were computed 
by first selecting host halos with V^°^ > 500 km s^^ in or- 
der to ensure the completeness of the subhalo count down 
to a threshold of V^^^^ > O.lbV^^lf (considering satellites 
above the scaled threshold roughly eliminates the scaling 
of satellite halo number with host halo size). We then 
fit the host density profile to an NFW profile [Eq. ^] 
in order to determine the best-fit Cvir- Finally, we used 
the subhalo counts and best-fit concentrations to compute 
(TVsat) and (cvir). The simulated host halos also show a 
clear correlation between concentration and satellite num- 
ber. The small number of host halos in this sample makes 
it difficult to draw detailed conclusions, yet the Cvir-Asat 



correlation measured in the simulation appears to be in 
good agreement with the model predictions. This corre- 
lation also appears in agreement with the related corre- 
lation between halo T^nax and s atellite number observed 
in the simulations of iGao et al.l (|2004al ). Fitting the sim- 
ulated halo results to a power law iVgat oc c^j" yields a 
slightly steeper relation than the model fit, with a best-fit 
power-law index of a « 0.98 ± 0.15. Given the errors, the 
power-law indices from the model and simulation data fits 
are consistent with each other. Interestingly, we find no 
statistically significant trend in the radial distribution of 
satellites with host halo concentration or collapse epoch. 

Figure[T2]shows the correlation between Oc and the num- 
ber of satellite halos selected according to a specific selec- 
tion threshold, namely satellites with V^mtx/KlJax ^ 0-15. 
However, we expect the strength of the correlation to be 
dependent upon the satellite sample threshold for several 
reasons. The drag of dynamical friction is proportional 
to Afg^ati so dynamical friction is most efficient at driving 
rather massive satellites toward the centers of their hosts. 
Furthermore, the median lookback time to a merger with a 
satellite of a specific mass is a decreasing function of Mgat- 
Mergers with large subhalos occur relatively more recently, 
on average, than mergers with small subhalos. In addi- 
tion, the accretion of large subhalos dominates the rate 
of host mass growth, so the accretion times of relatively 
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lookback' dyn lookback' dyn 

Fig. 11. — The competition between accretion and destruction at various redshifts. Left Panels: The thick solid lines show the accretion 
time distributions of subhalos with ViJ^' > 80kms~^ at the epoch of accretion as a function of lookback time to accretion for a host of mass 

1 q ^ _ 1 max — r- 

A^host = 10 h Mq. The lookback time has been scaled by the typical dynamical time of the host halo at the time of accretion, tdyn- The 
thin solid lines show the number of accreted subhalos that remain above the aforementioned Vmax threshold at the redshift of observation. 
From top to bottom, we show the model results at final redshifts of z; = 0, 1, 3. In all panels, the vertical dotted line depicts the median value 
of iio 

okback/^dyn foi" accreted subhalos. Right Panels: The fraction of surviving subhalos as a function of the lookback time to accretion. 
The lines show the fraction of halos that are above the threshold V^^^ ^ 80kms~^ at the time of accretion and remain above this threshold 
at the redshift of observation. Again, the lookback time has been scaled by the dynamical time of the host halo at the time of accretion. All 
panels show model results for a host of mass Mj^jst = 10^^'^ /i~^Mq at final redshifts of (from top to bottom) z = 0,1,3. Listed in each 
panel is the total number of subhalos that are accreted with V^^x ^ 80 km that remain above this threshold at the redshift of observation, 
^surv. The results shown in both panels are from 1000 realizations of the model of §[3] 
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Fig. 12. — Left Panel: Correlation between the collapse epoch as defined in Eq. ((Sj and the number of surviving subhalos. We plot the 
mean relation between the halo collapse epoch ac, scaled by the mean collapse epoch (flc), and the number of satellites Afgat, scaled by the 
mean number of satellites (Afgat), predicted by 1000 realizations of our substructure model. The central thick solid line corresponds to the 
mean relation for Mijo^t = 10^"* h~^MQ hosts. The other line types represent the mean relation at four different host masses: 10^^ h~^MQ 
{dot- dashed), lO" h-'^M^ (long- dashed), lO"-^ /i-^Mq {short- dashed), and lO" h-^MQ {dotted). The upper and lower thick, solid lines 
delimit the realization-to-realization scatter in the relation based on the 1000 model realizations of a Mhost = lO^'*'^ h~^'MQ host. We count 
only subhalos with VjJax ^ '^■^^Knax' sd, z = 0. Right Panel: The relationship between host halo concentration and the number of satellite 
halos. The lines represent the mean relation between halo Cvir scaled by the mean concentration (cvir), and the number of satellite halos 
scaled by the mean number of satellite halos (Afaat)i predicted by 1000 model realizations. The different line types correspond to different 
host masses and are the same as in the left panel. Again, only subhalos with V^ax ^ O-l^V^JJ are counted. The stars represent the measured 
concentrations and satellite numbers for all host halos in the ACDMgo simulation with V^^x* ^ 500kms~^. A power-law fit to the simulation 
points yields A^sat c~." with a ^ 0.98 ± 0.15. The relationship from the subhalo model yields a slightly smaller best-fit power-law index of 
a ^ 0.86 ± 0.07. 



massive subhalos more strongly influence any particular 
definition of formation epoch than the accretion times of 
small subhalos. These factors suggest that the strength of 
the flc-A^sat (or Cvw-Nsat) correlation should depend upon 
the size of the subhalos included in the sample and that 
the number of large subhalos should be more strongly de- 
pendent upon the mass accretion history of the host than 
the number of small subhalos. It is possible to use our 
substructure model to explore this threshold dependence, 
because our model provides a large number or realizations 
of host halos of various masses so that we are not ham- 
pered by the limited dynamic range of the cosmological 
simulation. 

We demonstrate the V^nax threshold dependence of the 
relation between iVgat and Uc in Figure 1131 The abscissa 
represents minimum values of scaled maximum circular 
velocity V^tx/^ax ' ^sed to define various subhalo sam- 
ples. At each threshold, we fit the number of satellite halos 
above the threshold to a power law 

Ns.t{> ^^ax/Kijax ) « «c- (10) 

At all thresholds, we find a power law to be an acceptable 
fit and plot the value of the best-fit power law indices on 
the vertical axis. The black squares show the relationship 
between b and (l^mtx/^max*) at z = 0. That the strength of 
the correlation between satellite number and Oc is a strong 
function of the T^iax threshold is evident in the Figure [T31 



The solid line shows the linear relation 

6 = 0.25-f4(y-ax/K'ax), (11) 

which we find to be a good fit to the relationship be- 
tween the power law index b, and the sample threshold 
K^ax/K^axS in the regime 0.08 < Kfax/K^ax < 0.5 at 
z = 0. As can be seen in Fig. [131 the strength of the 
A'sat-flc correlation declines with increasing redshift as dy- 
namical processes become less influential in determining 
subhalo populations. At z = 1, the best- fit relation for 
the power law index is 6 ~ '2.7{V,l^^/V,^lf). The marked 
evolution between redshift z = and z « 1 is due to the 
cosmological constant (recall 57a = 0.7) which causes a 
dramatic reduction in the rate of structure growth, and 
therefore satellite accretion rates, at low redshift. 

Another representation of the detailed dependence of 
the correlation between satellite abundance and host halo 
formation epoch (or host halo concentration) on the rela- 
tive subhalo size is displayed in Figure [Til for cluster-size 
host halos. The Figure shows the differential velocity func- 
tions (DVF) of subhalos as a function V^tx/^nax 
host halos and for two subsamples subdivided by host halo 
formation epoch. The curves in Fig. [TJl were constructed 
by stacking our 1000 model realizations of hosts at three 
masses, Mvh. = lO"-^, 10^^-^, and lO^-*-^ H'^Mq, for a to- 
tal of 3000 realizations, in order to overcome noise in the 
measurement. The solid line represents the mean DVF for 
all 3000 host halos. The upper, dashed line represents the 
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Fig. 13. — The Afsat-Qc relationship as a function of subhalo Vmax 
threshold. We show the best-fitting power law indices from a fit of 
the number of satellites above a particular Vmax threshold A^sat(> 
Knax/^max )' ^ power law function of the host halo collapse epoch 
Nss.t{> ^max/^ma?) ''c i ^ ^ function of the satellite threshold 
Knax/^max redshifts z = (squares), z = 1/2 (triangles), z = 1 
(circles), and 2 = 3 (stars). The solid line represents a fit of the 
power law index b, as a function of threshold to a line, b 2; 0.25 + 
4(K^ax/Knax*) at 2 = 0. The dashed line represents a linear fit at 

2 = 1, fe = 2.7(y-3^/yh-t). 

DVF for the half of the host halo sample with the highest 
Oc (lowest Cvir), while the lower, dot-dashed line repre- 
sents the DVF for the half of the sample with the lowest 
Uc- Again, the dependence of subhalo populations on host 
halo accretion history is evident. Late- forming, low con- 
centration host halos have more substructure of all sizes 
and larger subhalos are more strongly correlated with ac- 
cretion history and host halo Cvir for the aforementioned 
reasons. We speculate on possible implications of these 
relationships in § [SJ 

4.6. Application to Galaxy Clustering 

The model that we have presented has the virtues of 
being simple, quick and easy to compute, easy to parse and 
understand, and easy to modify or add specific ingredients 
to. These features make the model useful for studying 
a wide range of phenomena. In this section, we briefly 
discuss the particular example of applying this model to 
make predictions for the two-point correlation function of 
galaxies. This subject will be considered in much greater 
detail in the subsequent papers in subsequent papers of 
this series. 

We begin by using our model to compute the abun- 
dance of subhalos in a cosmological volume. By convolving 
the mean of the subhalo occupation distribution shown in 
Fig. El and Fig. [51 with th e known mass function of host 
halos () Jenkins et al.ll200l[ ). we can compute the number 
density of all halos, including both host halos and their 
subhalos, as a function of Knax threshold. The result is 




Fig. 14. — The dependence of the differential velocity function 
on host halo formation time and concentration. We show the differ- 
ential velocity functions of subhalos in cluster-size hosts computed 
by stacking the 1000 model realizations in three host mass bins, 
Afvir = lO^"' '*, 10"-5, and lO^*-^ H-'^Mq- The thick, solid central 
line represents the average differential velocity function for all 3000 
of the host halo realizations. The dashed line depicts the mean dif- 
ferential velocity function for the half of the host halo sample with 
the highest ac (lowest Cvir) and the lower, dot-dashed line repre- 
sents the differential velocity function for the half of the host halo 
sample with the lowest Oc (highest Cvir). In all cases, the error bars 
represent the estimated error on the mean value. 



shown in the bottom panel of Figure [TSl at three redshifts: 
z — 0, 1, and 3. The decrease in number density as a func- 
tion of redshift is a reflection of the fact that fewer massive 
halos have collapsed at earlier epochs. The dot-dashed line 
in the lower panel of Fig . 1 1 51 corresponds to the mean num- 
ber density of satellite halos only at z = 0. In the upper 
panel of Fig. [151 we show the fraction of halos that are sub- 
halos as a function of Vmax threshold, /sat(> Vmax)- At low 
redshift, /sat is roughly ^ 16 — 22% at thresholds that cor- 
respond to galaxy-size halos, which is comparable to the 
fraction of galaxies observed to be in groups and clusters. 
The decrease in /^at with increasing Vmax threshold reflects 
the fact that the mass function of host halos dn/dAfhost, 
is an increasingly steeply declining function of host halo 
mass. Therefore, as the threshold is increased, the num- 
ber density of host halos becomes increasingly dominated 
by host halos near the Vmax threshold. Halos are very un- 
likely to host satellites with a comparable Vnax (see Fig. [2]), 
so /sat decreases monotonically with increasing Vmax- A 
qualitatively similar argument applies to the dependence 
of /sat on redshift at fixed Vmax- The mass functions and 
the fraction of satellites are in good agreeme nt with results 
of simulations presented bv lKravtsov et ahl ([2004). 

In § [1] we suggested that our model of substructure 
could be coupled with a low-resolution, large - volume sim- 
ulation or with an analytic halo model (e.g. ISeliai3l200Clf ) 
in order to study the clustering properties of galaxies. We 
pursue these applications further in the forthcoming pa- 
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Fig. 15. — Halo and subhalo populations at three redshifts. Bot- 
tom Panel: The total number of halos (host halos and subhalos) 
with a maximum circular velocity greater than Vmax as a function 
of Knax. The solid line shows the halo population at 2 = 0, the 
dashed line shows the halo population at 2 = 1 and the dotted line 
shows the halo population at 2 = 3. The dot- dashed line shows the 
contribution from subhalos only at 2 = 0. Top Panel: The frac- 
tional contribution of subhalos to the total number of halos above 
each Vmax threshold at redshifts 2 = 0, 1, 3. The linotypes corre- 
spond to the same redshifts as in the bottom panel. 



pers of this series. However, the fact that the subhalo pop- 
ulations computed from our model are in good agreement 
with the subhalo populations of halos observed in cosmo- 
logical A^-body simulations does not imply that statistics 
that describe the spatial clustering of halos, like the two- 
point correlation function ^(r), will be in good agreement. 

Consider first, populating the host halos of a low- 
resolution simulation with subhalo populations derived 
from the model of § [3l In the most common imple- 
men tation of the EPS forma l ism for hierarchical cluster- 
ing (|Bond et all 119911 : ILC93I : ISomcrville & Kolatt 199§), 
which we have used to generate mass accretion histories 
for our host halos, the accretion history of a halo at fixed 
mass is independent of its large-scale environment. There 
is some evidence suggesting that this pr actical app roxi- 
mation may be justified (Lemson fc Kauffmannlll999t but 
see also iSheth fc Tor^^Ien[l2Q0l 7 Yet, significant correla- 
tions between environment and host halo properties, such 
as Oc or, more generally, the entire halo mass accretion his- 
tory, would also imply correlations between environment 
and subhalo populations. These correlations would then 
result in systematic differences between correlation func- 
tions measured directly from a simulated halo sample and 
those measured from a sample constructed by populating 



simulated hosts with subhalos from our model and it is 
important to understand this potential systematic effect 
before applying this method. As a byproduct, we also test 
the extent to which environmental effects may influence 
the correlation function and one of the fundamental as- 
sumptions of the halo model, the statistical independence 
of halo properties from halo environment. 

We attempted to estimate the potential importance of 
these effects using the following procedure. We identified 
host halos in the ACDMgo simulation as described in § [51 
and computed their virial masses and the positions of their 
most bound particles. We then used these masses and po- 
sitions as inputs to our substructure model. We computed 
a semi-analytic mass accretion history and a corresponding 
subhalo population for each host in the simulation, thereby 
constructing a hybrid catalog of simulated host halos and 
semi-analytic model subhalos. We reiterate that we did 
not use the concentrations or Knax values of the simulated 
hosts in order to construct the mass accretion histories 
and subhalo populations. Instead, we assigned each host 
halo new values of Cvir, Knax, and based on the semi- 
analytic mass accretion histories that we generated and 
the method outlined in § 13.21 As such, these host proper- 
ties and the properties of their subhalos are re-assigned in 
a way that is independent of any effects due to the large- 
scale environment. However, we did use the positions of 
the hosts to construct our hybrid catalogs, so the hybrid 
catalogs exhibit the same large-scale structures as the sim- 
ulation catalogs and there should thus be no differences 
between the two due to cosmic variance. Within each sim- 
ulated host halo, we place our model subhalos at a distance 
R/Rvir from the halo's most bound particle, where i?/i?vir 
is given by the dynamical model. The subhalos are thus 
forced to have spherical symmetry within their hosts and 
do not trace the radial profiles of the N-body halos they oc- 
cupy. We repeated this procedure ten times for each host 
halo in the ACDMgo simulation, thereby generating ten 
realizations of independent hybrid halo catalogs in order 
to account for variation from realization-to-realization of 
the semi-analytic model. Finally, we computed two-point 
correlation functions for the halos-|-subhalos in the hybrid 
catalogs and the ACDMgo simulation. 

We summarize the results of this experiment in Fig- 
ure ll61 where we show the correlation functions for samples 
selected above a Vmax threshold of Knax > 150kms~^. As 
a reference, this threshold corresponds to a mean number 
density of n « 8 x 10~^/i^Mpc~^ (see Fig. [TS]), which is 
close to the observed number density of galaxies with r- 
band absolute magnitudes of Af^^ — 19.5 ([Blanton et al.l 
I2003D . The squares in Figure [TBI represent the correlation 
function of the halos in the ACDMgo simulation. The error 
bars are the maximum of the jackknife error computed by 
excluding, in turn, each of the eight octants of the simula- 
tion volume and the Poisson error based on the number of 
pairs at each separation. A power-law fit to these points 
yields C(r-) ~ (r/5 h-^Mpc)-'^-^^ . 

The shaded band shown in Fig. [TBI shows the enve- 
lope of the ten correlations functions computed from our 
ten hybrid halo catalogs with host halo masses and po- 
sitions taken from the ACDMgo simulation and all other 
halo properties, particularly the subhalo populations, de- 
termined using the semi-analytic model. The correlation 
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functions of the halos of our hybrid model calculations are 
generally in remarkably good agreement with the corre- 
lation function of the simulated halos. On intermediate 
separations, r ^ 0.5 — 1 /i~^Mpc, ^(r) of the hybrid model 
halos is slightly higher 15%) than the correlation func- 
tion of the ACDMgo halos. These separations correspond 
to a scale where most halo pairs are from halos that re- 
side within the same host halo (including the host itself 
if it is part of the sample). In fact, the number of halo 
pairs on these scales is dominated by pairs within a com- 
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mon host in the mass range 10^"^^ -Mhost/ h ~^Mf. 
(for example, see iBerlind fc Weinberg! [2002f) . The num- 
ber of pairs within a single host halo is N{N — l)/2, 
where N is the number of halos within the host and 
iV = 1 + A^sat if the host is included in the sample while 
N = Nsat otherwise. Schematically then, ^(r) on separa- 
tions r ~ 0.5 - 1 h-^Mpc is roughly set by {N{N - 1)) 
for hosts with 10"< A/host/ h-^MQ< 10^'*. The slight en- 
hancement of ^(r) on these scales for the hybrid model is 
consistent with the differences in the second moment of 
the subhalo occupancy between the model results and the 
ACDMgo subhalos shown in Fig. [SI 

The fact that the two-point function of the hybrid model 
halos is in good agreement with the two-point function of 
the simulated ACDMgo halos is encouraging. First, the 
properties of our hybrid model halos are independent of en- 
vironment. Therefore, this calculation provides an explicit 
demonstration that any correlations of host halo proper- 
ties, such as Cvir, Oc, Nsat{> Knax); ^hc Spatial distri- 
bution of subhalos within their hosts, with environment 
are sufficiently weak as to be nearly unmeasurable in ^(r). 
This implies that populating comparably low-resolution, 
large-volume simulations with subhalos from our model is 
a viable method for studying the two-point statistics of 
galaxies, though one may likely use a more sophisticated 
mapping of galaxies onto halos and subhalos than the sim- 
ple Vmax threshold that we used to test the importance of 
environmental effects. We explore such models in a forth- 
coming paper. 

Additionally, t he hal o model of galaxy clu s tering 
(iNevman fc ScottI Il952t IScherrer fc Bertschinge^ Il99ll; 



ISeliakll2000l : IScoccimarro et al .Il2001h has been widely used 
recently for a variety of applications, one of which is to 
model galaxy clustering properties and to infer the galaxy 
occu pation of dark mat t er ha l os frorn obse rvational data 
fe.g. IZehavi et all 12001 [200i : IZheng||2004D . One of the 
fundamental assumptions of the halo model is that the 
galaxy occupation of dark matter halos is statistically in- 
dependent of host halo environment and depends only 
upon Mhost- Our results indicate that this is likely an 
acceptable assumption that does not lead to considerable 
systematic errors in £,{r). We close this section by illus- 
trating the utility of our model with an explicit, entirely 
analytic calculation of the two point function of halos us- 
ing the subhalo occupation distributions depicted in Fig. [5] 
and Fig. [51 coupled with the halo model. We compute 
the analytic c orrela tion function as described in detail in 
Zehavi et al. (|2004) using the host halo mass function of 



Jenkins et al.l ( 2001 ), the dark matter correlation function 



of [Smith et al.l (|2003( ). and the scale-dependent halo bias 
model of Tinker et al. (2004, in preparation). We ac- 
count for the radial profile of subhalos in our model by 
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Fig. 16. — Halo two-point correlation functions for all halos (host 
halos and subhalos) with Knax > 150kms~^. The squares represent 
the correlation function measured from the halos in the ACDMgo 
simulation. The error bars on the points represent the maximum of 
the jackknife error bars computed from the octants of the simulation 
box and the Poisson error based on the number of pairs. The shaded 
region represents the envelope spanned by ten realizations of our 
substructure model for each host halo in the simulation (see text 
ij l4.6l for details). The solid line represents the correlation function 
computed by using the halo occupation distributions in Fig. [5] and 
Fig. |6l in a purely analytic computation based on the halo model of 
clustering. 



fitting in each host mass bin, the mean radial distribu- 
tion of subhalos to a profile with a constant density core: 
nsati^) {1 + r/rc)~^. The resulting correlation function 
is shown as the solid line in Fig. [THj The agreement be- 
tween all three methods of computing ^(r) is apparent at 
all separations in the Fig. 1161 and is quite impressive. The 
fact that the halo model calculation is lower on large scales, 
8 h~^Mpc, is due to cosmic variance resulting from the 
finite size of the simulation box (L = 80 h^^Mpc). 

5. DISCUSSION AND CONCLUSIONS 

In this paper, we presented a semi- analytic model for 
subhalo populations within host dark matter halos. The 
details of the model arc discussed in § [31 with some com- 
plementary information in ZB03.. In the context of this 
model, the evolution of a subhalo population is treated 
as a multi-stage process and the transitions between each 
stage are treated as abrupt. The main ingredients of this 
model are as follows: 

1. Semi-analytic halo merger histories, computed using 
the EPS formalism; 

2. A method for assigning halo density profiles that in- 
corporates the known correlations between mass ac- 
cretion history and halo structure; 
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3. A distribution of initial orbital energies and angular 
momenta for merging subhalos; 

4. A prescription for modeling the orbital evolution of 
subhalos within hosts, including the effects of orbital 
decay via dynamical friction and tidal mass loss. 

In § [31 we make a detailed comparison between the re- 
sults of of our model and the results of cosmological TV- 
body simulations in the regime where they can be com- 
pared robustly. We calibrated our one-parameter model 
to match the normalization of the mean cumulative ve- 
locity function of subhalos within simulated, MW-size 
(Mhost ~ 10^^ /i-^Mg) host halos at z = 0, in a standard 
ACDM cosmology with a power spectrum normalization 
of (Tg = 0.9. Subsequently, we found, in the range that 
the simulation predictions are robust, our model to be in 
good agreement with the results of direct simulation. Wc 
summarize the success of our model when compared to 
simulation results with the following points. 

1. The velocity functions of cluster-size halos in a 
ACDM cosmology with 17a/ 1 - I^a = 0.3, h = 0.7, 
and fJs = 0.9 are in good agreement. 

2. We find that the mass functions of MW-size and clus- 
ter size halos in the same, ACDM cosmology are also 
in agreement. 

3. The model predicts a radial distribution of satellite 
halos within their hosts that is in good agreement 
with simulation results and exhibits a sizable anti- 
bias of subhalos with respect to the d ark matter. This 
is in co ntrast to the recent results of iTaylor fc Babul 
(|2004ch . 

4. We find a good agreement between the first and sec- 
ond moments of the halo occupation distribution pre- 
dicted by the model and measured directly in cosmo- 
logical simulations. The agreement holds for a wide 
range of host masses from Mhost ~ 5 x 10^^ /i~^Mq to 
Mhost ~ 3 X 10^'' H-^Mq at redshifts z = 0, 1 and 3 
in a ACDM cosmology with a normalization of erg = 
0.75, different from the normalization of the simu- 
lations used to provide inputs and to calibrate the 
model. 

5. We find a correlation between host halo collapse 
epoch Oc, and satellite number iVgat, at fixed host 
halo mass in our model. This implies a correspond- 
ing correlation between host Cvir and Ng^t at fixed 
Mvir. We find that A^^at oc c~-^, with a ~ 0.86 ± 0.07 
for subhalos selected such that T^mtx/KJSax* > 0-15. In 
the regime where this is testable with the simulations, 
we find a similar correlation. A fit to the simulation 
results yield a power law index a ~ 0.98 ± 0.15. Our 
model results suggest that this power law index varies 
rapidly with the subhalo selection criterion. 

6. Our model results in deviations from a self-similar 
scaling of subhalo abundance with host halo mass. 
In particular, we find that the number of satellite ha- 
los at a fixed value of Msat/Afhost scales in proportion 
to ~ -^^host' approxi mate agreement wit h the re- 
cent numerical study of 'Gao et al.' ('2004a') and the 
analytic work of Ivan den Bosch et al.. (,2004b ) . 



The success of our model in matching the results of sim- 
ulations in a variety of ways and at a variety of redshifts 
is non-trivial. It is also encouraging. This success sug- 
gests that such a model can be used for a wide range of 
applic ations. How e ver, n o te that in contrast to o ur re- 
sults, iTavlor et all (|2003D . iTavlor fc BabuH (|2004ct ). find 
that their analytic subhalo distributions are significantly 
more centrally concentrated than simulated subhalo popu- 
lations. The two models have several different ingredients 
and the root of this discord is not yet clear. Our results 
indicate that a centrally concentrated subhalo distribution 
is not a generic prediction of semi-analytic models. 

We also pointed out one possible shortcoming of the 
model. Although the semi-analytic procedure reproduces 
the mean quantities observed in simulations quite well, it 
seems that at least over some range of host properties, the 
model over-predicts the second moment of the subhalo oc- 
cupation (iVsat ( A'^jat — 1)) when (iVs at) is less than a few 
see Fig . [Hand Fig. O. The model of Ivan den Bosch et al.l 
2004bf) exhibits this same feature notwithstanding the 
fact that their treatment of mass loss is quite different 
from ours. We speculate that the origin of this discrep- 
ancy is a fundamental shortcoming in the standard EPS 
formalism for generating halo merger trees, althoug h fur- 
ther work is necessary to confirm this hypothesis. |W02| 
found that EPS merger histories have higher scatter in 
formation times than simulated merger histories (see also 
ILin. Jing. fcTh][2003h : given the relation found here be- 
tween iVgat and the formation time Oc, this discrepancy 
might have been anticipated. 

We demonstrated explicitly, in § 14. 4[ that deviations in 
the self-similar scaling of subhalo populations with host 
halo mass are due to a relative shift in the balance of 
power between accretion rates and the orbital evolution 
timescales of subhalos. This argument has a natural ex- 
tension. In ? I4.5[ we took advantage of the fact that semi- 
analytic approaches are not subject to the same restric- 
tions on sample size as direct simulation, to expound upon 
correlations between host halo properties and their satel- 
lite populations. 

We found that the host halo collapse time Oc correlates 
strongly with the number of satellite halos. The sense 
is such that early-forming halos have fewer satellites sim- 
ply because they acquired their satellites earlier and these 
satellites were therefore subject to the destructive process 
of the dense, host environment for a longer time. This is 
in qualitative agre e ment w ith the recent, similar analysis 
of ITavlor fc Babull (l2004bD . who use a different definition 
of halo formation time. Format ion time is strongly cor - 
related with halo concentration (jW02l : IZhao et"ani2003l ). 
implying that satellite halo abundance is correlated with 
halo concentration. In our model, host halo concentra- 
tion varies in inverse proportion to Oc via Eq. |4l so this 
correlation occurs by construction, with A^gat c";^'*^ for 
sateUites selected by V^fjV^l^ > 0.15. 

We confirmed that this correlation between halo con- 
centration and subhalo abundance is exhibited by host ha- 
los in simulations (Fig. [T^ . For the ACDMgo simulation, 

we find that A^,at(l^rax/'Cax > 0.15) cx c;i°■9^ in good 
agreement with the results of our model. An extension of 
this argument is that the strength of the flc-iVsat (or Cvir- 
A'sat) correlation should depend upon the Knax threshold 
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used to define the satellite halo sample. We find that at 
any particular threshold of ^max/Knax ' -^sat oc aj, where 
b 0.25 + 4(T/-yFhost)^ < y^^jV^l^ < 0.5 
at redshift z = 0. We find that the correlation between 
formation time and satellite abundance (or, equivalently, 
between host concentration and satellite abundance) be- 
comes weaker at high redshift as shown in Fig. [131 By 
redshift z = 1, the relation between the number of satel- 
lites and the host halo formation time is described by a 
power law index 5 ~ 2.7(y-yKSoBt). 

These results may have direct observational implica- 
tions. For example, cluster halos of fixed mass but with be- 
low average substructure counts are objects that acquired 
their substructure earlier, allowing satellites to experience 
significant dynamical evolution. If accretion into a cluster 
environment strips gas from the subhalo and effectively 
truncates star formation, this suggests that these early- 
forming clusters may have a higher fraction of red galax- 
ies. More directly, in a sample of clusters restricted about 
a narrow range of X-ray temperatures or luminosities, clus- 
ter richness may be expected to correlate with the fraction 
of red galaxies. 

Similarly, central galaxy mergers would be more com- 
mon in these early-forming systems, likely causing the 
central galaxy to become more luminous. This may in- 
duce a correlation between cluster richness and the lu- 
minosity of the brightest cluster galaxy in systems of 
fixed dynamical mass. We are tempted to identify the 
so-called fossil groups with the extreme, early-formation 
tail of the formation epoch (ac) distribution in our model. 
The fossil groups are systems that are characterized by X- 
ray luminosities comparable to those of poor galaxy clus- 
ters, while their optical luminositie s are dominated by a 
single, bright, cen tral galaxy (e.g., IVikhlinin et al.|[l999t 
I Jones et"al1 12003|) . In this scenario, these fossil groups 
would correspond to large, group-size halos that assem- 
bled their mass very early, when dynamical timescales were 
short, allowing any substructure to evolve significantly. 
Any large satellites would then be more vulnerable to a 
rapid merger with the central object due to dynamical 
friction, leaving a luminous central galaxy, and to severe 
tidal disruption, leaving a surrounding halo bereft of lu- 
minous companions. This interpretation is consistent with 
the analyses of the radial distributio n of low-lu minosity 
satellite galaxies in such groups (Mathews et al. |[2004.) . 

The correlation between A'^gat and Cvir may be detectable 
observationally. Our results suggest that for clusters with 
similar mass estimates, the optically-richer clusters should 
have underlying dark matter halos with concentrations 
that are lower than average. It may be possible to detect 
directly trends of this kind using cluster mass profile esti- 
mates either via X-ray analyses or gravitational lensing. If 
a trend between formation time or concentration and the 
number of satellite galaxies Ngai , could be detected using 
optical data, it could potentially provide an optical mass 
measurement with less scatter than Ngai or total cluster 
luminosity. Specifically, a large fraction of the scatter in 
mass at a fixed observed value of Ngai in the newest gen- 
eration of optically-selected cluster catalogs is due to the 
theoretically-expected scatter in Ngai at fixed mass (R. H. 
Wechsler et al., in preparation). The aforementioned cor- 
relation suggests that clusters at fixed Ngai are a mix of 



later-forming, low-mass clusters and earlier-forming, high- 
mass clusters. Thus, it may be possible to remove some of 
this scatter with other optical observables. 

Finally, we demonstrated the utility of the presented 
model with an explicit calculation of correlation functions 
in § 14.61 We generated a set of ten hybrid halo catalogs 
using the host halos from the ACDMgo simulation and ten 
semi-analytically computed subhalo populations for each 
host. We computed the two-point correlation function of 
the ACDMgo halos and compared it with the correlation 
functions of the halos in the hybrid catalogs. We found 
the correlation functions calculated in this way to be in 
good agreement. The host concentrations and the hybrid 
model subhalos are assigned according to a technique that 
is ignorant of the environment of the halo. This result 
implies that any correlations between Cvir, A"sat(> Knax)i 
or subhalo spatial distributions with environment are suf- 
ficiently weak that they likely do not have a measurable 
effect on the correlation function. This supports one of the 
basic as sumptions of halo model a nalysis of galaxy cluster- 
ing fe.g. lZehavi et al.ll2003l I2004D . namely that the galaxy 
occupation of host halos is statistically independent of en- 
vironment and depends solely upon the mass of the halo. 
This also suggests that the technique of populating a large- 
volume, low-resolution simulation with semi-analytic sub- 
halo populations is viable for studying galaxy clustering 
and environment. We pursue these issues in a forthcom- 
ing paper in which we also study mapping galaxies onto 
halos using various prescriptions in order to test several 
hypotheses about the physics of galaxy formation. 
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